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INTRODUCTION 


There  are  many  classes  of  boundary  value  problems  which  are  not  effectively  solved  using  a  single 
numerical  method.  These  problems  might  be  more  effectively  solved  by  combining  different 
methods,  applying  each  method  to  portions  of  the  domain  where  it  is  best  suited.  In  this  study  we 
consider  the  coupling  of  the  boundary  element  method  (BEM)  with  the  finite  element  method 
(FEM)  to  solve  problems  in  elastostatics. 

Boundary  element  methods  have  become  an  accepted  alternative  to  domain-based  methods,  such  as 
finite  element  and  finite  difference  methods,  for  many  classes  of  bout  dary  value  problems  (Refs  1 
and  2)  The  direct  boundary  elment  method  (DBEM)  and  indirect  boundary  element  method 
(IBEM)  are  two  well  established  formulations  that  are  used  to  solve  a  variety  of  boundary  value 
problems  in  engineering.  The  key  to  these  formulations  is  the  existence  of  a  fundamental  singular 
•'olution  (i.e.,  the  solution  for  a  point  source  or  point  load  in  an  infinite  domain  -  a  free-space 
Green's  function). 

The  fundamental  soludon  is  jointly  the  blessing  and  curse  of  the  methods.  Since  it  is  a  solution  of 
the  governing  differential  equations,  problems  defined  by  a  self-adjoint  differential  operator  can  be 
expressed  by  integral  equations.  With  the  governing  equations  satisfied  only  the  boundary 
conditions  remain  to  be  satisfied.  Numerically  this  means  that  for  many  classes  of  problems  BEMs 
require  only  the  boundary  to  be  discretized.  The  dimensionality  of  these  problems  is  reduced  by 
one.  Tlie  methods  thus  lend  themselves  to  modeling  problems  with  infinite  domains  since  only  the 
surface  of  cavities,  not  the  domain,  must  be  di.scretized.  This  decrease  in  dimensionality  reduces 
the  number  of  degree^  of  freedom  (dof)  in  the  approximation;  however,  the  resulting  equations  are 
unsymmetric  and  fully-populated  in  contrast  to  the  symmetric  sparse  systems  characteristic  of 
domain-based  methods.  The  singular  nature  of  the  fundamental  solution  is  suited  to  solving 
problems  having  high  stress  gradients.  Yet  this  same  singular  nature  complicates  the  numerical 
integrations  inherent  to  the  methods  when  higher  order  elements  are  employed.  Fundamental 


singular  solutions  exist  for  many  classes  of  boundary  value  problems  with  homogeneous  domains. 
However,  for  an  arbitrary-inhomogeneous  domain  the  fundamental  solution  is  not  known,  and 
additional  numerical  effort  is  lequired  (Ref  3).  Since  the  methods  are  relatively  new,  they  have  not 
been  as  extensively  developed  as  domain-based  methods  for  nonlinear  applications.  They  also  lack 
the  generality,  in  terms  of  extensive  continuum  and  structural  element  libraries,  that  commercial 
finite  element  computer  programs  possess. 

The  finite  element  method  (FEM)  has  been  the  most  highly  developed  approximate  solution 
technique  for  many  classes  of  boundary  value  problems  (Ref  4).  The  most  general  derivation 
follows  from  a  weighted  residual  approximation  of  the  boundary  value  problem.  The  error 
introduced  by  an  assumed  approximate  solution  form  is  orthogonalized  with  respect  to  a  set  of 
weighting  functions  to  obtain  the  "best"  approximate  solution.  For  the  Galerkin  weighted  residual 
statement  the  weighting  functions  correspond  to  the  basis  functions  in  which  the  approximate 
solution  is  expanded.  For  problems  with  a  self-adjoint  operator,  the  weak  form  of  the  Galerkin 
statement  is  equivalent  to  an  alternative  variational  statement  of  the  problem. 

Basis  functions  in  the  FEM  usually  are  locally  based  polynomials  for  which  numerical  integrations 
are  relatively  simple.  The  local  nature  of  these  functions  leads  to  a  sparse  system  of  equations,  but 
this  same  attribute  hinders  its  ability  to  tnodel  infinite  domains.  To  compensate  for  this  limitation, 
Ungless  (Ref  5,  1973)  and  Bettess  (Ref  6, 1977)  developed  special  elements  (infinite  elements)  to 
model  infinite  domains  which  include  shape  functions  "based  on"  the  form  of  analytical  solutions 
(i.e.  they  should  decay  radially  outward  towards  infinity  and  integrations  over  the  domain  should 
remain  finite).  The  generality  of  the  Galerkin  and  variational  statements  ahows  nonlinearities  and 
domain  inhomogencities  to  be  accommodated.  For  most  structural  applications  the  resulting 
system  of  equations  is  also  symmetric. 

The  boundary  and  finite  element  methods  are  related.  In  seeking  to  model  problems  having  infinite 
domains,  Zienkiewicz,  Kelly  and  Bettes  (Ref  7,  1977)  interpreted  the  IBEM  as  a  generalized  finite 


element  methodology  where  the  "shape  functions"*  were  obtained  directly  in  terms  of  the 
fundamental  solution.  The  DBEM  can  be  interpreted  in  the  same  manner  but  the  unknown 
coefficients  are  nodal  boundary  tractions  and  displacements  instead  of  artificial  boundary  tractions. 
(The  physical  identification  of  the  coefficients  made  here  is  only  valid  for  "continuous  elements" 
using  collocation.)  This  contrasts  with  Brebbia's  woric  (Ref  8,  1978)  relating  the  analytical  basis 
of  the  methods  (FEM,  DBEM,  and  finite  difference  method)  where  he  derives  Somigliana's  first 
identity  by  employing  the  fundamental  singular  solution  as  the  weighting  function  in  the  weak 
statement  of  the  boundary  value  problem. 

Considering  the  individual  strengths  of  boundary-based  and  domain-based  methods,  some  classes 
of  problems  may  be  more  effectively  solved  by  combining  the  methods.  Nonlinear  soil-structure 
interaction  and  fracture  mechanics  are  two  classes  of  problems  where  a  combined  solution 
approach  can  potentially  be  more  effective.  For  nonlinear  soil-structure  interaction  problems,  the 
FEM  could  be  used  to  model  the  structure  and  a  region  of  the  soil  which  is  expected  to  exhibit 
nonlinear  behavior,  and  the  BEM  could  be  used  to  approximate  the  infinite  domain.  For  brittle 
fracture  mechanics  problems,  the  FEM  could  be  used  to  model  all  of  the  problem  except  a  region 
local  to  the  crack  tip  where  the  BEM  could  provide  a  crack  tip  "element"  capable  of  representing 
very  high  stress  gradients.  (While  BEMs  have  been  formulated  for  modeling  inelastic  material 
behavior,  the  coupling  approach  of  this  study  assumes  the  BEM  region  is  linear  elastic.  Ductile 
fracture  could  be  solved  by  a  coupled  approach  but  maybe  with  no  gain  in  effectiveness.)  For  both 


*  The  term  "shape  function"  has  been  used  in  various  ways  by  researchers  working  with  coupled  solution 
approaches.  Many  uses  of  the  term  in  referring  to  the  BEM  are  a  bit  of  a  misnomer.  In  the  initial  part  of  the  work 
by  Zienkiewicz  et  al.  (Ref  7)  "shape  function"  is  defined  as  a  function  which  spans  the  BEM  region  (a  "super 
element"),  but  unlike  typical  FEM  shape  functions  it  is  in  general  nonzero  at  all  nodes.  In  this  context  the  shape 
functions  provide  a  basis  for  the  displacement  field  however  the  variation  of  the  shape  function  along  a  given  portion 
of  the  boundary  is  not  known  a  priori  to  be,  for  example,  linear  or  quadratic.  For  the  BEM,  "shape  function"  is 
usually  used  to  refer  to  a  boundary  interpolant  used  along  a  so  called  "boundary  element"  in  performing  the  numerical 
integrations.  These  "shape  functions"  are  nonzero  at  a  single  node;  however  they  do  not  form  a  basis  for  the 
boundary  response,  only  for  the  boundary  distribution  assumed  for  the  integrations.  This  is  discussed  in  more  detail 
in  a  later  section.  This  later  definition  will  be  used  in  this  study. 


of  these  problems  one  aspect  of  the  problem  is  infinite:  for  the  soil  problem  the  domain,  for  the 
crack  problem  the  stresses.  While  one  might  argue  that  both  of  these  are  mathematical  fictions, 
they  do  effectively  characterize  the  behavior  of  the  respective  problems. 

Objective 

The  objectives  of  this  thesis  are  to:  (1)  formulate  sti^ness  matrices  for  both  the  indirect  and  direct 
boundary  element  methods  that  would  allow  them  to  be  coupled  to  the  finite  element  method,  (2) 
develop  a  coupling  algorithm  which  requires  a  minimum  of  modification  to  existing  separate  FEM 
and  BEM  software  systems,  (3)  implement  one  of  the  stiffness  matrix  formulations,  and  (4)  study 
the  basic  behavior  of  the  coupled  solutions  including  accuracy  and  the  incompatibility  along  the 
FEM-BEM  interface.  The  ultimate  objective  of  this  research  is  to  determine  whether  a  combination 
of  boundary  element  and  finite  element  solution  methods  could  more  effectively  solve  nonlinear 
structural/geotechnical  problems. 

Scope 

In  this  study,  I  am  investigating  different  approaches  to  obtaining  a  FEM-hosted  coupling.  The 
emphasis  is  on:  (1)  conceptually  unifying  some  of  the  different  coupling  techniques,  (2)  providing 
an  algorithm  to  facilitate  implementation,  and  (3)  studying  the  numerical  behavior  of  a  coupled 
system.  A  physically  intuitive  development  of  an  IBEM  stiffness  matrix  is  included.  This  direct 
derivation  of  a  stiffness  matrix  is  followed  by  the  corresponding  variational  derivation.  Though 
the  emphasis  is  on  the  IBEM,  the  DBEM  is  also  addressed  The  numerical  study  is  limited  to  two 
classical  problems  which  illustrate  the  behavior  of  the  coupled  system.  Overviews  of  both  the 
DBEM  and  IBEM  are  given  in  following  sections  to  establish  nomenclature  and  to  present  the 
necessary  equations  for  background. 
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Background  on  Coupled  Solution  Approaches 

The  idea  of  coupling  the  FEM  and  BEM  is  attributed  to  Wexlcr  (Ref  9, 1972),  who  used  integral 
equation  solutions  to  constrain  FEM  solutions  for  infinite  domain  field  problems  in  the  early 
1970s.  Otlicr  early  examples  of  using  a  combined  solution  approach  dealt  with  wave  phenomena, 
Chen  and  Mei  (Ref  10,  1974)  and  Shaw  (Ref  11,  1974).  Zienkiewicz  et  ai  (Ref  7,  1977) 
proposed  a  combined  solution  approach  for  statics  problems  followed  by  Osias,  Wilson  and 
Seitelman's  (Ref  12,  1977)  u.se  of  a  coupled  solution  to  solve  elastostatic  problems. 

There  are  numerous  approaches  to  coupling  the  methods  (Ref  13).  In  this  study  the  main 
cla.ssification  is  ba.sed  upon  "which  method  hosts  the  coupling."  For  a  BEM-hosted  coupling  the 
FEM  subdotnain  is  treated  as  a  BEM  region;  equilibrium  and  compatibility  are  approximately 
enforced  continuously  along  the  interface,  analogous  to  the  BEM  approtich  to  modeling  piece- wise 
homogeneous  problems.  The  resulting  equations  arc  unsymmetric  and  thus  this  approach  only 
lends  itself  to  problems  numerically  dominated  by  the  BEM  subdomain  (Refs  13,  14,  15,  16  and 
1 7).  For  a  Fl'.M  hosted  coupling  the  BEM  subdomain  is  treated  as  one  or  many  finite  elements.  In 
this  "super-element"  approach  the  resulting  BEM  equations  reflect  the  "equivalent"  stiffness  of  the 
sulxlomain  and  can  lie  directly  assembled  by  a  FEM  program.  The  term  "super-element"  reflects 
the  high  degree  of  connectivity  which  will  typically  be  present  for  a  BEM  region.  A  coupled 
solution  approach  does  not  seem  to  necessitate  a  displacement  based  FEM;  however  due  to  its 
commercial  success,  it  appears  all  the  studies  of  coupled  solutions  have  been  based  on  this 
Ibrmulation. 

In  addition  to  the  HEM  hosted  and  FEM-hosted  coupling  approaches,  there  have  been  a  number  of 
approac  hos  wlterc  in  a  sense  neither  methixl  hosts  (as  defined  above)  the  coupling.  This  open- 
ended  caiei’i'iy  mainly  ineludes  approaches  where  one  methotl  serves  as  a  boundary  condition  for 
the  other,  .Shaw's  study  (Ref  1 1,  1974)  of  time  hiurnonic  scattering  is  an  early  example  of  such  an 
approaeh. 


Shaw  and  Falby  (Ref  18, 1977)  devised  another  approach  in  which  the  boundary  integral  equation 
was  incorporated  into  a  variational  formulation  for  the  FEM  region  as  a  natural  boundary  condition 
(also  see  Ref  19). 

The  treatment  of  the  BEM  as  a  boundary  condition  to  the  FEM  region  has  also  been  accomplished 
iteratively  (Ref  20).  The  BEM  regions  are  initially  treated  as  if  they  are  rigid  and  the  FEM  problem 
is  solved.  The  reaction  forces  along  the  interfaces  of  the  regions  are  applied  to  the  BEM  regions 
and  the  corresponding  boundary  value  problems  for  the  BEM  regions  are  solved.  The 
displacements,  obtained  by  the  BEM  analyses,  along  the  FEM-BEM  interface  are  applied  to  the 
FEM  region  and  the  problem  is  re-solved.  After  a  few  iterations  the  interfaces  converge  to  their 
equilibrium  position. 

Kamat  and  Brown  (Ref  21,  1985)  proposed  yet  another  iterative  approach.  Instead  of  using  a 
BEM,  an  "ancestor"  to  the  IBEM  is  used  where  the  unknowns  are  the  positions  and  strengths  of  a 
finite  number  of  point  sources  (for  the  potential  problem).  Their  approach  is  "scalar"  in  nature  in 
the  sense  that  it  does  not  require  the  explicit  formulation  of  a  stiffness  matrix.  The  potential  energy 
functional  for  the  FEM  region  is  augmented,  via  a  penalty  number,  by  a  compatibility  condition 
along  the  FEM-BEM  interface.  This  compatibility  condition  requires  a  discrete  least  squares 
satisfaction  (i.e.,  at  three  times  as  many  points  along  the  interface  as  there  are  point  sources).  The 
augmenferl  functional  is  then  minimized  using  a  preconditioned  conjugate  gradient  technique. 

Within  FEM-hosted  couplings  the  first  level  of  classification  is  the  procedure  used  to  obtain  the 
stiffness  relationship.  The  two  approaches  are:  (1)  a  variational  approach  where  a  boundary 
variational  principle  is  the  basis  of  the  relationship,  and  (2)  a  direct  approach  where  the  BEM 
equations  are  manipulated  in'o  a  stiffness  form.  A  key  step  to  both  approaches,  using  either  BEM 
formulation,  is  establishing  a  relationship  between  the  nodal  tractions  and  the  nodal  displacements; 
an  exception  to  this  is  a  variational  approach  with  the  IBEM  which  retains  the  artificial  tractions  as 
unknowns  in  the  system  of  equations  (Refs  7, 13  and  14). 


Chen  and  Mei  (Ref  10,  1974)  were  apparently  the  first  to  formulate  a  coupled  solution  based  on  a 
variational  principle.  In  their  study  of  the  action  of  sea  waves  on  a  harbor  their  "BEM"  region  was 
expressed  as  a  finite  series  of  Bessel  functions.  They  used  a  modified  functional  for  the  FEM 
region  to  enforce  continuity  of  the  dependent  variables  and  their  derivatives  across  the  FEM-BEM 
interface  as  a  natural  boundary  condition.  The  outer  "BEM"  region  was  referred  to  as  a  "super 
clement."  Other  early  efforts  using  a  variational  approach  include:  Jeng  and  Wexler  (Ref  22, 
1977):  Zienkiewicz,  Kelly,  and  Bettess  (Ref  7,  1977);  and  Kelly,  Mustoe,  and  Zienkiewicz  (Ref 
13,  1979). 

Approximation  of  the  variational  principle  necessarily  produces  a  symmetric  system,  and  thus  these 
approaches  are  sometimes  grouped  with  some  direct  approaches  as  "symmetric  couplings."  In  the 
variational  approach  considered  in  this  study,  a  variational  statement  of  the  problem  is  written  in  a 
boundary  form  with  the  absence  of  body  sources  or  forces  (Refs  7  and  23).  For  elastostatics  the 
functional  corresponds  to  a  boundary  form  of  the  potential  energy  functional.  There  have  been 
many  modifications  of  this  functional  for  various  purposes  including  (1)  to  relax  essential 
boundary  conditions  or  interface  compatibility  and  (2)  to  enforce  equilibrium  (Refs  23  and  24). 

In  approximating  the  variational  statement,  relationships  between  boundary  values  obtained  from 
the  integral  equations  are  substituted  into  the  functional,  reducing  the  number  of  unknowns.  For  a 
dispF'.cement-based  FEM,  the  nodal  tractions  arc  usually  expressed  in  terms  of  the  nodal 
displacements,  which  requires  the  inversion  of  a  fully-populated  coefficient  matrix.  The 
previously  mentioned  exception  to  this  approach  (using  the  IBEM),  writes  the  boundary 
displacements  and  tractions  in  terms  of  the  artificial  tractions.  In  this  case  stationarity  is  taken  with 
respect  to  both  the  nodal  displacements  (of  the  adjacent  "framing"  region)  and  the  nodal  artificial 


boundiiry  tractions  instead  of  nodal  displacements  alone.  The  lack  of  a  matrix  inversion  is  paid  for 
by  an  increase  in  the  number  of  unknowns*  and  a  double  integration  process. 

In  contrast  to  the  variational  formulations,  Brebbia  (see  e.g.,  Ref  2,  8,  15  or  25)  directly 
manipulated  the  BEM  equations  into  a  stiffness  form.  In  general  direct  formulations  do  net 
produce  a  symmetric  stiffness  matrix.  Brebbia  symmetrized  the  system  by  averaging  the  off- 
diagonal  terms  based  on  an  error  minimization  argument  using  the  method  of  least  squares.  The 
resulting  stiffness  matrix  was  of  the  same  form  as  obtained  by  variational  methods  using  the 
potential  energy  functional  (Ref  7). 

Despite  the  apparent  agreement  between  the  two  approaches,  agreement  between  proponents  of  the 
approaches  was  scarce.  Kelly  et  al.  (Ref  13)  stated  that  the  modification  to  obtain  symmetry 
"apparently  has  no  valid  foundation  except  that  the  energy  approach  led  to  a  similar  form." 
Brebbia  el  al.  (Ref  2)  later  concluded  that  there  was  no  reason  for  the  stiffness  matrix  to  be 
symmetric  other  than  for  reasons  of  computational  convenience  and  efficiency.  In  addition,  he 
questioned  the  combination  of  variational  principles  and  integra  l  equations  used  in  the  variational 
approach.  Initial  effon.s  indicated  that  the  symmetric  matrix  provided  better  results;  however,  more 
recent  comparisons  (Ref  24,  1982)  using  quadratic  elements  and  providing  for  geometric 
discontinuities  indicate  the  contrary.  Fortunately  for  cases  which  have  a  relatively  fine  boundary 
descretization,  Roudolphi  (Ref  26)  has  found  the  degree  of  asymmetry  to  be  insignificant.  This  is 
likely  to  be  the  case  for  geotechnical  applications. 

Hartman  (Ref  27)  gives  a  somewhat  humorous  overview  of  the  problem  as: 

"lliis  boundary  element  method,  as  it  is  called,  leads  to  nonsymmetric  and  fully- 
populated  matrices  and  it  is  not  ea.sy  to  convince  engineers,  that  anything  good  can 
come  out  of  such  strange  matrices.  Imagine  these  same  engineers  when  they 


Zicnkicwicz  et  al.  (Rd  7)  eliminate  the  coefficients  corresponding  to  the  artificial  boundary  traction  disUibution  at 
the  clement  Icvc'.  Zicnkicwicz  (Refs  4  and  13)  also  points  out  the  direct  analogy  that  exist  bctw'een  this  particular 
formulation  and  that  of  hybrid  finite  elements. 


discover  that  the  stiffness  matrices  which  they  derive  from  these  linear  systems  turn 
out  to  be  nonsymmetric... .  God  knows  there  is  not  much  engineers  believe  in,  but 
they  certainly  believe  that  stiffness  matrices  have  to  be  symmetric  and  they  are 
therefore  deeply  irritated,  finding  themselves  in  a  situation  where  two  fundamental 
beliefs  clash  (Mathematics  and  Mechanics)  and  either  one  is  equally  convincing." 

Another  problem  which  is  common  to  both  approaches  is  failure  of  the  stiffness  matrix  to  satisfy 
equilibrium.  Equilibrium  requires  the  terms  in  each  column  corresponding  to  a  given  global 
direction  to  sum  to  zero.  Boissenot,  Lachat,  and  Watson  (Ref  28,  1974)  augmented  the  system  of 
equations  for  the  DBEM  with  two  additional  equations  which  cause  the  approximate  traction 
distribution  to  identically  satisfy  equilibrium  (in  English  see  e.g.  Refs  13  and  23).  By  this 
approach  the  equilibrium  problem  is  addressed  while  establishing  the  nodal  traction-displacement 
relationship  instead  of  after  the  stiffness  matrix  is  formed.  Hartman  (Ref  27)  gives  a  theoretical 
discussion  on  the  symmetry  and  equilibrium  problems  and  provides  advice  on  how  to  correct 
them.  Tullberg  and  Bolteus  (Ref  24)  compare  the  accuracy  and  convergence  characteristics  of 
seven  stiffness  formulations  for  two  simple  two-dimensional  problems  —  uniaxial  tension  and 
bending.  The  seven  stiffness  matrices  were  obtained  using  the  DBEM  and  included  all  the 
variations  mentioned  above.  In  their  study  the  direct  formulation,  in  it’s  unsymmetric  form,  gave 
the  best  accuracy  and  rates  of  convergence. 

In  addition  to  the  differences  used  here  to  classify  formulations,  there  have  been  numerous  other 
variations.  Some  researchers  (Ref  23)  have  used  element  condensation  to  reduce  the  unknowns 
for  the  BEM  super-element  to  those  located  on  the  interface  between  the  two  representations.  This 
could  be  useful  in  fracture  mechanics  applications  where  numerous  boundary  elements  might  be 
used  near  the  crack  tip.  As  previously  mentioned,  modified  functionals  have  been  used  to  relax 
compatibility  requirements  (e.g.,  when  interfacing  subdomains  with  different  "shape  functions"). 
As  will  be  shown  in  this  thesis,  even  when  the  "shape  functions"  between  methods  match  the 
interfaces  are  inherently  incompatible. 

Variations  in  the  BEM  formulations  themselves  introduce  additional  classification  parameters. 
Though  not  addressed  in  this  study,  for  some  applications  traction  discontinuity  is  an  important 
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consideration.  Many  ways  of  treaing  this  problem  have  been  reported  (Refs  24  and  29-31).  A 
simple  iliough  less  rigorous  approach  is  to  posiuon  the  collocation  points  at  the  interior  of  the 
element.  When  all  of  the  collocation  points  associated  with  an  element  are  internally  positioned, 
the  elements  are  referred  to  as  "discontinuous"  or  "nonconforming"  elements  (Ref  32).  Though 
the  "discontinuous"  elements  add  some  flexibility  to  the  method  in  terms  of  mesh  refinement,  it  is 
at  the  expense  of  increasing  the  number  of  equations  by  up  to  100  percent.  Discontinuous 
elements  would  appear  to  complicate  a  coupling  algorithm  since  the  collocation  points  do  not 
correspond  to  node  points. 

Relationship  to  Previous  Work  by  the  Author 

The  study  this  thesis  is  based  on  was  funded  by  the  Naval  Civil  Engineering  Laboratory  (NCEL). 
This  section  gives  an  overview  of  the  work  which  led  to  the  current  study.  In  our  research  at 
NCEL,  the  primary  role  of  the  BEM  has  been  to  compliment  FEM  analysis  capabilities  (i.e., 
coupled  solution  techniques).  I'he  BEM  is  in  general  a  more  analytically  complicated  technique, 
and  thus  the  fundamentals  of  the  method  required  investigation.  As  a  result  of  these  efforts  the 
potential  of  the  method  as  a  separate  analysis  tool  has  also  been  observed.  Thus  our  couplii.^ 
efforts  have  emphasized  a  modular  approach  where  it  is  assumed  that  both  analysis  methods  can 
serve  as  an  individual  or  a  coupled  tool. 

In  our  first  study  (Ref  33)  we  investigated  the  IBEM  formulation  in  one-  and  two-dimensional 
elastostatics  and  compared  numerical  results  with  the  DBEM  and  analytical  solutions.  We  used  a 
crude  element  (constant  distribution)  but  the  potential  of  the  method  was  evident.  This  study 
indicated  the  need  for  higher  order  elements  and  integration  techniques  to  improve  results  in  the 
near  boundary  region. 

In  our  next  study  (Ref  17)  we  investigated  a  coupling  of  the  BEM  and  FEM.  We  formulated  a 
BEM-hosted  coupling,  approximating  compatibility  and  equilibrium  on  the  interfaces. 
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Qualitatively  the  coupling  was  a  success,  but  the  simple  constant  elements  used  for  both  methods 
prevented  substantial  quantitative  evaluation.  For  nonlinear  soil-structure  applications  where  the 
FEM  is  used  to  model  the  nonlinear  region,  the  FEM  should  host  the  coupling.  Two  areas 
required  further  study:  (1)  accuracy  of  the  BEM  formulations,  and  (2)  a  FEM-hosted  coupling 
approach. 

In  our  third  study  (Refs  34  and  35)  we  investigated  higher  order  isoparametric  boundary  elements 
and  their  associated  integration  techniques.  Lachat  and  Watson  (Ref  36)  adapted  the  isoparametric 
element  formulation  directly  from  FEM  technology.  However,  for  the  BEM  the  integrands  in  the 
coefficient  calculations  are  singular  and  thus  their  integration  requires  special  attention.  Though 
much  effort  had  been  devoted  to  de.aling  with  the  integration  over  the  singularity  an  effective 
numerical  approach  to  performing  the  integrations  in  the  near  boundary  region  appeared  to  be 
lacking.  We  developed  a  recursive  algorithm  which  adaptively  subdivides  the  element  to  maintain 
consistent  accuracy  independent  of  the  response  point  location.  The  new  element  formulation  has 
effectively  eliminated  near  boundary  error  resulting  from  coarse  numerical  integration  and  has 
provided  valuable  insight  to  other  errors  inherent  to  the  boundary  element  methods.  This 
technique*  is  the  key  to  allowing  a  numerical  demonstration  of  the  incompatibility  which  occurs  in 
the  coupling  (see  the  numerical  results  section). 

The  initial  effort  in  this  study  was  presented  in  a  preliminary  report  (Ref  38).  This  thesis 
completely  supersedes  the  previous  report. 


*  This  technique  serves  as  a  "shell"  to  standard  numerical  integration  techniques  like  Gauss  quadrature.  A  dramatic 
increase  in  efficiency  is  attainable  by  improving  the  quadrature  rule  directly.  A  "shell"  might  still  be  necessary  but  it 
would  be  resorted  to  less  frequently.  For  a  recent  paper  on  improved  integration  methods  see  reference  37. 
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Nomenclature 

This  is  a  list  of  principal  symbols  used  within  the  thesis.  Additional  variables  which  have  a  limited 
scope  of  use  are  defined  within  the  text.  Modifications  (such  as  the  addition  of  a  subscript)  are 
also  defined  in  the  text  unless  given  below.  Instead  of  the  matrix  notation  defined  below,  indicial 
notation  will  sometimes  be  used  to  clarify  equations.  Indices  will  always  be  lower  case  while 
identifying  subscripts  or  superscripts  will  ’  e  upper  case. 

The  most  general  form  of  a  variable  (i.e.,  with  the  least  amount  of  supplementary  symbols)  is 
defined.  For  example,  u^  represents  a  displacement  vector  but  depending  on  the  notation  it  may  be 
continuous  or  discrete,  known  or  unknown.  There  are  some  non-unique  uses  of  Latin  symbols. 
For  these  cases  supplementary  symbols  are  included  as  necessary.  Indicial  notation  is  only  used  in 
the  definitions  below  for  tensor  quantities;  the  index  range  in  this  case  is  two.  In  the  following 
definitions  "V"  represents  a  generic  variable. 

Mathematical  Symbols 

V(x)  Explicitly  defines  the  quantity  as  continuous  as  opposed  to  discrete. 

V  A  rectangular  or  square  matrix  of  functions  (e.g.,  M(x))  or  constants  (e.g.,  K). 

V  A  column  vector. 

A 

V  A  known  value(s);  for  example,  u  is  a  vector  of  known  nodal  displacements. 

V  An  unknown  value(s);  for  example,  l  is  a  vector  of  unknown  nodal  tractions. 
Transpose  of  a  matrix  or  vector.  is  a  row  vector. 

v-i 


Inverse  of  a  matrix. 
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Latin  Symbols 


A^^  A  constant  tensor  which  provides  a  zero  displacement  for  the  fundamental 

singular  solution  at  a  reference  distance  r^.  Ajj^  =  -5iicCj[C2ln(ro)  -1], 


C  Square,  sparse  matrix  which  contains  boundary  integrations  of  traction  and 

displacement  "shape  function"  products. 

Cl,  C2,  C3,  C4  Constants  used  in  defining  the  fundamental  solution.  Functions  of  the  material 
constants  (see  text). 


e)c(D  ^  unit  concentrated  force  used  in  the  definition  of  the  fundamental  solution. 

E  Young's  modulus. 


E 

Eijk(i.4) 


f 

F 


Gij(x,D 

G 


MOO 


N(20 


n 


Square  matrix  relating  nodal  tractions  to  nodal  displacements,!  =  E  u. 

Kernel  function  in  Somigliana's  2”*^  identity  relating  boundary  displacements  to 
the  stress  field. 

Column  vector  of  generalized  nodal  forces. 

Kelvin  solution  for  traction. 

A  matrix  of  element  integrations  involving  the  Kelvin  solution  for  traction. 
Both  the  direct  and  indirect  formulations  consist  of  F  matrices.  The  two  F 
matrices  differ  due  to  the  switch  in  roles  of  the  field  and  source  point  and  the 
indices. 

Fundamental  singular  solution,  Kelvin  solution  for  displacement 

A  matrix  of  element  integrations  involving  the  fundamental  solution.  Due  to  the 
symmetries  in  the  fundamental  solution  the  G  matrices  for  the  direct  and  indirect 
formulations  are  equivalent. 

Kernel  function  in  Somigliana’s  2"**  identity  relating  boundary  tractions  and 
body  forces  to  the  stress  field. 

Traction  "shape  functions"  used  with  nodal  traction  values  to  approximate  a 
piecewise  continuous  traction  distribution. 

Displacement  "shape  functions"  used  with  nodal  displacement  values  to 
approximate  a  continuous  displacement  distribution. 

Order  of  the  system  of  equations. 
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Pj  Artificial  boundary  traction  vector,  analogous  to  simple-layer  source  or  source 

density  of  potential  problems. 

r  Distance  from  the  field  to  the  source  point,  r^  =  y^j. 

tj  Actual  boundary  traction  vector. 

Tijk(x,4)  Kelvin  solution  for  stress. 

Uj  Displacement  vector. 

x;  Position  vector  to  a  point  on  the  boundary  or  within  the  domain  of  the  problem 

depending  upon  the  context. 

yj  Vector  from  the  source  to  the  field  point,  y,  =  x;  - 

Zi  Position  vector  for  a  point  in  the  domain  used  in  the  integration  of  body  forces. 

Greek  Symbols 

r  Complete,  finite  boundary  of  the  problem. 

5ij  Kronecker  delta  symbol, 

ejj  Cartesian  strain  tensor. 

X  A  Lame  constant. 

|i  A  Lame  constant,  the  shear  modulus. 

V  Poisson's  ratio. 

Position  vector  to  a  point  on  the  boundary  or  within  the  domain  of  the  problem 
depending  upon  the  context. 

ajj  Cartesian  stress  tensor. 

Q  Domain  of  the  problem. 

¥i 


Vector  of  body  forces  per  unit  volume. 
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THEORY 


This  study  focuses  on  the  development  of  a  stiffness  matrix  for  a  coupled  solution  approach.  Each 
method  of  analysis  is  applied  to  portions  of  the  domain  where  it  is  best  suited.  For  illustration  the 
domain  Q  (Figure  1)  is  divided  into  two  subdomains  d**,  the  FEM  subdomain,  and  d®,  the  BEM 
subdomain.  The  equations  of  equilibrium  are  given  by 


< 

II 

O 

(la) 

a.  =  o.. 

U  J» 

(lb) 

and  the  strain-displacement  relations  are  given  by 


e..  =  —  (u.  +  u. .) 

ij  2 


(2) 


The  BEM  subdomain  is  assumed  to  be  a  linear  isotopic  homogeneous  elastic  material.  Thus  the 
governing  constitutive  relations,  generalized  Hooke's  law,  are  given  as 

<j,j  =  2>iey  +  (3) 


where  X  and  p.  are  the  Lame  constants  expressed  in  terms  of  Young's  modulus  (E)  and  Poisson's 
ratio  (v)  as 


X  = 


Ev 


(l-hv)(l-2v) 


Vi  = 


2(1 +v) 


(4) 


Alternatively  the  equilibrium  (1),  strain-displacement  (2),  and  constitutive  (3)  relations  can  be 
combined  to  give  a  single  set  of  equations  in  terms  of  displacement  written  as 

uu. +  (X  -I-  u.)u. ..  -I-  w.  =0 


(5) 


'6 


the  Navier  equations. 

The  boundary  conditions  are  given  by 

U;(x)  a  (20 

xer^j 

(6a) 

II 

a  £  r 

(6b) 

\vb'  '  ^ii(20  and  ti(x)  are  presoribed  distributions  of  boundary  displacements  and  tractions, 
respectively.  The  simple  notation  does  not  imply  that  the  boundary  conditions  are  mutually 
exclusive;  the  fully-mixed  boundary  value  problem  is  addressed. 

The  key  solution  to  equation  (5)  for  the  more  popular  BEMs  is  the  fundamental  singular  solution  -- 
obtained  for  elasticity  by  Kelvin.  This  is  the  solution  due  to  a  concentrated  load  in  an  infinite  space 
which  has  the  same  dimension  as  the  problem  to  be  solved.  For  plane  strain  conditions  the  Kelvin 
solution  gives  the  responses  due  to  a  unit  force  e^C^)  in  an  infinite  plane.  The  indices  i,  j,  and  k 
assume  values  of  1  or  2,  and  repeated  indices  imply  summation.  The  Kelvin  solution  gives  the 
displacement  field  as 

UjCx)  =  Gj,^(x,^e^(^  (7a) 

where 


87tp.(l  -  v) 

C2  =  3  -  4v 

Ajk  =  arbitrary  constant  tensor  based  on  zero  displacement  reference  distance 

Yi 


17 


=  ViVi 

By  incorporating  equations  (2)  and  (3),  the  stress  field  o^OO  is  given  as 

Oj.OO  =  T..^(2i^e^(^  (8a) 

where 


C4  =  1  -  2v 

Equilibrium  conditions  applied  at  a  boundary  point,  with  a  unit  outward  normal  njOO,  and  equation 
(8a)  combine  to  give  the  surface  tractions  tjCs)  as 

ti(20  = 


where 


Figures  2  and  3  illustrate  the  singular  behavior  of  the  Kelvin  solution  components  Gn  and  Tm 
respectively,  for  a  point  load  applied  at  the  origin  of  the  Cartesian  system.  The  singularity  of  G  is 
order  ln(r),  while  the  singularity  of  T,  obtained  from  derivatives  of  G,  is  order  1/r.  The  plane 
strain  solution  can  be  convened  to  the  plane  stress  solution  by  specifying  an  effective  Poisson  ratio 
v  =  v/(l+v) 
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The  fundamental  solution  is  a  key  ingredient  in  formulating  integral  equations  for  the  indirect  and 
direct  BEMs.  Integral  equations  are  an  equivalent  statement  of  a  boundary  value  problem.  Finite 
difference  solutions  approximate  the  differential  equations;  finite  element  solutions  approximate  a 
variational  principle  or  weak  statement  of  the  boundary  value  problem;  and  BEMs  approximate  the 
integral  equations. 

The  following  subsections  provide  an  overview  of  both  the  direct  and  indirect  BEMs  and  their 
associated  integral  equations.  These  overviews  are  followed  by  derivations,  both  direct  and 
variational,  of  an  IBEM  stiffness  matrix.  The  last  subsection  discusses  the  discontinuous  nature  of 
coupled  solution  techniques  -  an  inherent  trait  of  integral  equation  methods. 

Overview  of  Direct  Boundary  Element  Method 

The  direct  boundary  element  method  is  the  most  highly  developed  of  all  the  integral  equation 
methods.  It  was  first  applied  to  elastostatics  by  kizzo  (Ref  39).  Cruse  and  Rizzo  (Ref  40) 
followed  with  a  solution  of  the  general  transient  elastodynamic  problem.  Since  these  early 
applications,  it  has  matured  numerically  and  the  scope  of  application  has  significantly  broadened. 

The  integral  equations  on  which  the  DBEM  is  based  are  known  as  the  Somigliana  identities.  The 
first  identity  has  been  obtained  by  two  different  derivations.  In  both  cases  the  fundamental 
solution  plays  a  key  role.  The  first  derivation  uses  Betti's  reciprocal  work  theorem  in  which  one 
system  is  the  actual  boundary  value  problem  and  the  second  system  corresponds  to  the 
fundamental  solution  (Ref  39).  The  second  derivation,  due  to  Brebbia  (Ref  8,  1978),  starts  with 
the  weak  statement  of  the  boundary  value  problem  and  incorporates  the  fundamental  solution  as  the 
test  or  weighting  function.  Betti  (1872-73)  and  Somigliana  (1885-86)  were  the  first  to  apply 
potential  methods  to  elasticity  (Ref  23);  thus  it  is  not  surprising  that  Betti's  theorem  and 
Soniigliana’s  first  identity  are  the  Navier  equation  counterparts  to  Green's  second  and  third 
fomiula  from  potential  theory. 
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Somigliana's  first  identity  is  given  by 

up  =  l|(is)  -  Fij&i)  U,(i)]  ds  +  K(ii)  >{>1 00  dj  (10a) 

r  n 

(a  e  r  A  2  e  n) 

in  which  the  displacement  field  is  written  in  terms  of  the  boundary  tractions,  boundary 
displacements,  and  body  forces.  The  second  identity  is  obtained  by  combining  equation  (10a)  with 
the  strain-displacement  (2)  and  constitutive  (3)  relations  giving  the  stress  tensor  as 

t  (JS)  -  Uj(20  j  dA  +  V.  (2)  d2  (10b) 

r  Q 

in  which  the  stress  field  is  now  expressed  in  terms  of  boundary  tractions,  boundary  displacements, 
and  body  forces.  Since  the  strain-displacement  relations  (2)  involve  derivatives  of  displacement, 
the  kernel  functions  H  and  E  have  singularities  of  orders  1/r  and  lA^,  respectively. 

The  DEEM  is  formulated  by  the  numerical  approximation  of  equations  (10).  Equation  (10a)  is 
used  to  obtain  the  unknown  boundary  values,  and  then  both  equations  allow  the  calculation  of 
internal  responses.  The  two  approximations  made  to  obtain  the  unknown  boundary  values  are:  (1) 
integrating  in  a  piecewise  manner  and  (2)  solving  the  equation  in  a  boundary  collocation  sense. 
That  is,  the  integration  is  subdivided  over  boundary  elements  and  domain  cells;  and  the  integral 
equations  are  applied  to  a  discrete  number  of  points  on  the  boundar>'. 

The  early  applications  of  the  method  used  analytical  integrations  over  elements  which  could  model 
constant  or  linear  distributions  of  boundary  and  domain  values.  More  recent  applications  apply 
isoparametric  element  concepts  common  to  the  finite  element  method  (Ref  36).  However,  unlike 
the  FEM  all  quantities  are  interpolated  at  the  same  order  (i.e.,  linear  boundary  elements  interpolate 
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geometry,  displacements,  and  tractions  lineariy).  The  boundary  distribudon  of  any  field  variable  y 
on  a  single  element  is  written  as 


X  N  (n) 

a  KX 


(11) 


where  a  is  the  node  index;  N£t4  is  the  number  of  element  nodes;  T)  is  the  normalized  local 
curvilinear  coordinate;  Yio  values  of  the  field  variable;  and  Nq(t))  are  the  appropriate 

polynomial  "shape  funcdons"  for  the  element.  A  typical  quadradc  element  =  3)  is  shown  in 
Figure  4.  The  "shape  functions"  for  this  element  are  given  by 


N,(t1)  =  2 

(12a) 

N2(ti)  =  (1-ti^ 

(12b) 

N3(ri)  =  j  (Tl^  +  Tl) 

(12c) 

The  kernel  functions  in  the  integral  equations  are  singular  and  require  special  treatment  when 
considering  the  response  at  points  on  (Ref  1)  or  near  (Ref  34)  the  boundary.  The  corresponding 
numerical  details  are  omitted  here  for  brevity.  While  the  use  of  the  fundamental  solution  ensures 
the  satisfaction  of  the  governing  differential  equations,  the  boundary  collocation  ensures  violation 
of  the  boundary  conditions.  This  point  is  discussed  in  a  later  section  with  respect  to  a  coupled 
solution  approach. 

The  piecewise  integration  and  boundary  collocation  allow  equation  (10a)  to  be  approximated  by  a 
system  of  linear  equations  in  terms  of  the  nodal  boundary  values  and  body  forces  as 

Gl  -  Fu  =  Q 


(13) 
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(for  details  see  Ref  1  or  2).  The  body  forces  are  specified  and  thus  the  product  with  G'  gives  u 
known  vector.  For  a  mixed  boundary  value  problem  the  vectors  of  nodal  tractions  and 
displacements  contain  known  and  unknown  values.  To  solve  for  the  unknown  nodal  values  we 
partition  the  system  as 


[G.  G,]j 

f.l 

k] 

1  - 

w 

1  ■* 

w 

+  G‘]ir  -Q 


(14) 


and  then  collect  the  unknown  nodal  values  to  give 


(15) 


The  unknown  nodal  values  are  obtained  by  solving  the  above  system.  Then  Somigliana's 
identities  can  be  used  to  obtain  the  desired  internal  responses. 


Overview  of  Indirect  Boundary  Element  Method 

The  indirect  boundary  element  method  is  probably  the  second  most  common  integral  equation 
method.  Like  the  direct  method  it  also  has  its  origins  in  classical  potential  theory  (Refs  4 1 , 42  and 
43).  Single-  and  dou’oa  layer  potentials  were  used  in  the  theory  of  classical  electrodynamics  to 
express  boundary  value  problems  as  integral  equations.  The  most  advanced  implementation  of  the 
method  appears  to  remain  in  its  application  to  electromagnetic  field  problems  (Refs  1 1, 23  and  44). 
The  indirect  method  appears  to  have  first  been  applied  to  elastostatics  by  Massonnet  (Ref  45, 
1965).  The  numerical  development  of  IBEM,  for  elasticity  problems,  paralleled  that  of  the  DBEM 
until  the  mid  197()s.  In  recent  years  it  has  not  been  developed  as  extensively  as  the  direct  methoil; 
however,  its  physically  meaningful  formulation  provides  insight  to  b  ihods  (Ref  35). 
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The  integral  equation  on  which  the  IBEM  is  based  can  be  derived  from  Somigliana's  first  identity 
by  considering  two  displacement  boundary  value  problems  associated  with  the  domain,  Q,  and  its 
compliment  space  (i.e.  the  interior  and  exterior  problems  with  respect  to  F).  Baneijee  and 
Butterfield  (Ref  1)  present  the  analogous  development  for  steady  state  and  transient  potential  flow 
based  on  an  idea  originally  due  to  Lamb  (Ref  46).  Eringen  and  Suhubi  (Ref  47)  derive  the  integral 
equations  for  elastodynamics  by  the  same  argument.  The  brief  presentation  of  the  IBEM  which 
follows  is  based  on  a  physically  intuitive  argument. 

The  integral  equations  on  which  the  IBEM  is  based  are  a  single-layer  potential  statement  of  the 
boundary  value  problem.  The  domain  is  embedded  in  an  infinite  plane  as  shown  in  Figure  5. 
For  elasticity  applications  the  single-layer  source  corresponds  to  a  vector  of  artificial  tractions, 
Pk(&)  ^  e  r.  In  this  formulation  we  seek  the  boundary  distribution  of  artificial  tractions  which 

satisfy  the  prescribed  boundary  conditions.  These  tractions  are  artificial  in  the  sense  that  they  only 
exist  because  the  domain,  Q,  has  been  embedded  in  an  infinite  plane.  They  represent  an 
intermediate  step  in  the  formulation  and  do  not  correspond  to  the  actual  boundary  tractions. 
Artificial  tractions  and  body  forces  can  be  expressed  as  a  "continuous"  distribution  of  e|j(^).  Thus 
the  displacement,  stress  and  traction  fields  are  given  by  integrating  equations  (7a),  (8a)  and  (9a) 
respectively  as 

Ui(x)  =  JGj^(x,^)  P^(4)  d^  +  |Gj^(x,z)  ijl5^(z)  dz  (16) 

r  ti 

o,j(x)  =  jTijk(x,^  Jt..^(x,z)  \(z)  dz  (17) 

r  £j 

t(x)  =  jF,k(2i,l)  Pk(^)  d^  -h  jFy^(i,z)4(i)dz 
r  n 


(18) 
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Since  the  solution  is  expressed  as  a  superposition  of  the  fundamental  solution,  the  governing 
differential  equations  are  satisfied  over  the  entire  plane  (excluding  F). 

To  determine  the  distribution  of  the  artificial  tractions  we  bring  the  field  point  x  to  the  boundary  F 
and  enforce  the  boundary  conditions  (6).  The  resulting  integral  expressions  are  given  as 

Ui(x)  =  Jgj^(x,^)  d^  +  |Gj^(x,z)  ^5^(2)  dz  xeFu  (19) 

r  £1 

t(x)  =  ^  ^  (20) 

r  a 

Equation  (19)  is  regular  upon  integration  while  equation  (20)  must  be  interpreted  in  a  Cauchy 
principal  value  sense  and  is  thus  written  as 

t.(xd  =  ±^5ikPk(l0  + +  jFj^(x.2)^(z)dz  zeF^  (21) 

r  n 

The  first  term  represents  the  singularity  contribution  assuming  a  tangent  line  through  x;  the  sign 
depends  on  the  orientation  of  the  element  with  respect  to  £2.  The  boundary  integral  must  be 
interpreted  as  a  Cauchy  principal  value  integral. 

The  IBEM  is  formulated  by  the  numerical  approximation  of  equations  (16)  through  (21). 
Equations  (19)  and  (21)  are  used  to  determine  the  unknown  artificial  tractions,  and  then  equations 
(16)  through  (18)  allow  the  calculation  of  internal  responses.  As  with  the  DBEM,  the  two 
approximations  made  to  obtain  the  unknowns  are:  (1)  integrating  in  a  piecewise  manner  and  (2) 
solving  the  integral  equations  in  a  weighted  residual  sense  (on  the  boundary).  The  integration  is 
subdivided  over  boundary  elements  and  domain  cells.  Normally  the  integral  equations  are  satisfied 
in  a  collocation  sense;  however.  Lean  et  al.  (Ref  44)  report  improved  accuracy  by  a  Galerkin 
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approximation."  This  study  is  limited  to  a  collocation  approximation  of  the  equations.  Equations 
(19)  and  (21)  are  then  approximated  as 

li=GE  +  G'\jl:  (22a) 


1  =  FE  +  F'4 


(22b) 


where  ^  and  ^  are  values  of  boundary  displacements  and  tractions  at  collocation  points.  For 
"continuous  elements"  the  collocation  points  correspond  to  element  node  points.  The  coefficient 
matrices  G,  G',  F,  and  F'  are  obtained  by  integrations  of  (19)  and  (21)  with  respect  to  the 
appropriate  "shape  functions."  £  is  the  vector  of  unknown  nodal  artificial  traction  values,  and  the 
last  term  provides  the  effect  of  the  body  forces  at  each  collocation  point.  For  details  see  reference 
1,  2,  35  or  54. 


For  a  mixed  boundary  value  problem^  equations  (22a)  and  (22b)  are  applied  to  Fy  and  r-]- 
respectively  giving 


A 


M 


(23) 


The  nodal  artificial  tractions  are  obtained  by  solving  the  above  system.  Then  equations  (16) 
through  (18)  can  be  used  to  obtain  the  desired  internal  responses.  As  when  developing  the  system 
of  equations,  boundary  values  are  obtained  by  letting  a  go  to  the  boundary  (as  (19)  and  (21)). 


*  For  additional  information  on  the  application  of  the  Galerkin  and  related  approximations  to  integral  equations  see 
Refs  7,  48.  49,  50,  51  and  52.  Over-collocation  has  also  been  applied  in  the  approximation  of  integri  equations 
(Ref  53). 

+  Improved  IBEM  formulations  have  been  developed  to  address  the  mixed  boundary  value  problem  (Ref  55)  which 
incorporate  dislocation  doublets  as  well  as  artificial  tractions  as  unknowns. 
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Consider  a  few  of  the  differences  between  the  two  BEMs.  Let  N^p  be  the  number  of  boundary 
collocation  points.  Ignoring  the  integrations  associated  with  body  forces,  the  DBEM  (15)  and 
IBEM  (23)  require  the  calculation  of  2(2Ncp)^  and  (2Ncp)^  coefficients,  respectively.  However, 
in  solving  the  system  of  equations,  the  DBEM  yields  the  unknown  boundary  values  directly  while 
the  IBEM  yields  the  artificial  traction  distribution.  To  calculate  all  the  unknown  boundary  values 
the  IBEM  requires  an  additional  (2Ncp)^  coefficient  calculations.  Though  the  total  number  of 
coefficients,  and  thus  integrations,  required  to  obtain  the  unknown  boundary  values  are  equal  --  the 
computational  effort  is  not.  The  coefficient  matrices  of  the  DBEM  are  calculated  "simultaneously" 
and  thus  the  "overhead  calculations"  (e.g.,  calculation  of  the  Jacobian)  associated  with  the 
numerical  integrations  are  done  once.  Since  these  "overhead  calculations"  are  performed  twice  by 
the  IBEM,  it  requires  more  effort.  However,  this  comparison  has  assumed  the  analyst  needs  all 
the  boundary  values.  They  are  required  by  the  DBEM  to  calculate  internal  responses  by  the 
Somigliana  identities;  they  are  not  needed  by  the  IBEM  for  internal  response  calculation  and  thus 
the  analyst  can  be  selective. 

The  calculation  of  internal  responses  also  differs  between  the  methods.  The  DBEM  obtains 
internal  responses  by  the  Somigliana  identities  (10)  which  integrate  the  effects  of  both  the 
boundary  tractions  and  displacements.  The  IBEM  only  integrates  the  effects  of  artificial  tractions, 
(16)  through  (18).  As  in  developing  the  system  of  equations,  though  the  effort  is  not  doubled  for 
the  DBEM,  it  is  considerably  more.  Another  important  difference  in  the  internal  response 
calculations  is  the  order  of  singularity  of  the  kernel  functions.  For  two  dimensional  problems  the 
methods  have  the  following  orders  of  singularity: 

Displacement  Strain.  Stress.  Traction 

DBEM  o(lnr)ando(l/r)  o(l/r)  and  o(l/r^) 

o(lnr)  o(l/r) 


IBEM 
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Thus  the  DBEM  must  deal  with  stronger  singularities  and  correspondingly  more  difficult 
integrations.  This  problem  is  most  severt  he  near  boundary  region  (Refs  34  and  35).  In 
fairness,  the  indirect  method  is  not  without  its  problems  in  calculating  internal  responses.  Though 
the  integration  effort  is  reduced  in  the  IBEM,  often  the  accuracy  is  too.  A  problem  that  is  inherent 
to  the  IBEM  is  associated  with  loading  and  geometric  discontinuities.  In  these  areas,  the  artificial 
tractions  experience  very  hig..  ^  iients  even  though  boundary  conditions  may  be  very  well 
behaved.  Thus  unless  these  areas  receive  special  treatment  in  the  numerical  foratulation  or 
modeling,  accuracy  is  locally  very  poor. 

For  coupled  solution  approaches,  particularly  for  infinite  domain  problenv ,  .t  is  not  easy  to 
determine  which  BEM  would  be  the  most  effective.  Both  methods  must  calculate  two  coefficient 
matrices  which  involve  the  same  orders  of  singularity.  In  addition,  discontinuities  can  often  be 
avoided,  eli.ninating  them  as  a  factor.  In  the  following  section  the  IBEM  is  used  to  illustrate  a 
direct  derivation  of  a  stiffness  matrix.  The  ideas  can  be  applied  to  the  DBEM  but  are  explained 
with  the  IBEM  due  to  its  conceptual  simplicity. 

A  Direct  Derivation  of  an  IBEM  Stiffness  Matrix 

This  section  presents  a  physically  intuitive  derivation  of  an  IBEM  stiffness  matrix.  The  relation 
between  nodal  displacements  and  tractions  is  identical  to  the  relation  obtained  directly  by  Kelly  et 
al.  (Ref  13).  The  physical  approach  provides  additional  insight  to  the  stiffness  formulation. 

A  stiffness  matrix  relates  nodal  displacements  to  generalized  nodal  forces  in  the  form 

Kii  =  f  (24) 

where  K,  n  and  ?  are  the  stiffness  matrix,  displacement  vector,  and  generalized  force  vector, 
respectively.  This  basic  definition  and  a  numerical  solution  provided  by  the  IBEM  allow  a 
stiffness  matrix  to  be  calculated  for  a  BEM  region.  Consider  the  displacement  boundary  value 
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problem  shown  in  Figure  6.  The  boundary,  F,  is  subdivided  into  Ng  boundary  elements.  The 
elements  are  shown  as  isoparametric  quadratic  elements  for  illustration.  All  displacements  are  zero 
except  the  j***  dof  which  has  a  unit  displacement.  We  can  now  use  the  IBEM  to  solve  the 
displacement  boundary  value  problem.  As  given  by  equation  (22a)  with  body  forces  omitted,  the 
unknown  nodal  artificial  boundary  tractions  are  related  to  the  known  nodal  displacements  by 

U  =  GE  (25) 

This  system  of  equations  can  be  solved  for  the  nodal  artificial  traction  values  which  can  then  be 
used  to  obtain  the  unknown  real  boundary  tractions  by  an  approximation  of  equation  (18).  The 
nodal  tractions  can  be  expressed  in  matrix  form  as 

1  =  (26) 

The  distribution  of  the  traction  along  the  boundary  is  then  approximated  in  terms  of  the  locally 
based  "shape  functions"*  as 

l(2D  =  M(s)i  (27) 

which  by  incorporating  equation  (26)  becomes 

l(x)  =  M(20Fi  (28) 

With  the  distribution  of  boundary  tractions,  generalized  nodal  forces  (on  the  boundary)  can  now  be 
obtained.  The  nodal  forces  are  determined  by  "measuring"  the  work  done  by  the  tractions  due  to  a 
series  of  virtual  displacements  —  displacement  "shape  functions."  The  generalized  forces  are  given 
by 


*  In  this  context  the  "shape  functions"  are  locally  based  but  defined  globally.  That  is,  they  are  nonzero  on  one  or 
two  elements  but  defined  on  all  f. 


28 


£  =  Jn^Oi)  lOO  dx 
r 


(29) 


For  the  prescribed  displacements  u  the  generalized  forces  are  the  j*  column  of  a  stiffness  matrix 
for  the  domain  Cl.  To  completely  calculate  the  stiffness  matrix  we  must  individually  subject  each 
remaining  nodal  dof  to  a  unit  perturbation,  solve  the  corresponding  displacement  boundary  value 
problem,  and  then  determine  the  generalized  forces  by  equation  (29)  to  obtain  each  column  of  K. 


We  now  seek  to  express  the  above  procedure  in  equation  form.  Using  indicial  notation,  the  nodal 
displacement  vector  for  the  J*  perturbation  is  given  by 


1,  i  =  J 
0,  i#J 


(30) 


where  i  ranges  from  1  to  the  number  of  dof  (n).  Each  is  the  J*  coordinate  vector  in  the  n 
dimensional  space.  Combining  all  of  the  vectors  into  a  single  matrix  (i.e.,  the  superscript  J 
becomes  an  index  j)  allows  equation  (25)  to  be  rewritten  as 


u.. 

ij 


(31) 


By  equation  (30),  Uy  is  simply  the  identity  matrix.  Thus  we  can  easily  solve  for  the  artific'i'tl 
tractions  as 


u  y 


(32) 


This  formulation  provides  insight  to  the  character  of  G'^;  each  column  of  G  ®  is  the  vector  of 
artificial  tractions  resulting  from  a  unit  perturbation  of  the  corresponding  dof.  G  is  of  full  column 
rank  and  thus  invertible  if  the  collocation  point  positions  are  unique.  Inserting  this  result  into 
equations  (28)  and  (29)  yields  the  stiffness  matrix  as 
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K,j  =  Jn^iOO  1^,00  (33) 

r 

or  in  matrix  notation 

K  =  Jn’^OO  M(x)  dx  F  G'^  (34) 

r 

Kelly  et  al.  (Ref  13)  obtained  the  same  form  of  solution  for  the  potential  problem  by  eliminating  the 
source  density.  For  the  elastostatics  problem  this  is  equivalent  to  eliminating  £  from  equations 
(25)  and  (26).  With  their  approach  we  see  that  the  product  of  F  and  G"^  provides  the  relation 
between  nodal  tractions  and  displacements.  We  can  write  equation  (34)  in  an  abbreviated  form  as 

K  =  CE  (35) 

where 

C  =  Jn'^Cs)  MCs)  di  (36) 

r 

and 

E  =  FG''  (37) 

Brebbia's  direct  approach  (Ref  25)  to  obtaining  a  stiffness  matrix  for  the  DBEM  also  has  the  same 
form.  For  the  DBEM  however,  E  is  given  by 

E  =  G'^  F  (38) 

where  F  is  not  the  same  as  obtained  for  the  IBEM.  The  DBEM  uses  the  same  fundamental 
solution;  however  the  roles  of  the  source  and  field  points  and  the  roles  of  the  indices  are  reversed 
in  obtaining  the  Somigliana  identities. 
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As  previously  mentioned,  direct  formulations  of  stiffness  matrices  in  general  do  not  initially  yield  a 
symmetric  matrix.  The  following  variational  derivation  of  an  IBEM  stiffness  matrix  somewhat 
justifies  taking  the  symmetric  component  of  the  above  K. 

A  Variational  Derivation  of  an  IBEM  Stiffness  Matrix 

This  section  presents  a  variational  derivation  of  an  IBEM  stiffness  matrix.  The  boundary 
variational  statement  on  which  this  derivation  is  based  was  presented  by  Zienkiewicz  et  al.  (Ref  7). 

As  an  alternative  to  the  previously  given  strong  statement  of  the  boundary  value  problem  an 
equivalent  variational  statement  can  be  used.  The  solution  to  the  boundary  value  problem  is  the 
displacement  field  which  satisfies  the  essential  boundary  conditions  (6a)  and  minimizes  the 
potential  energy  functional  given  by 

n(u)  =  -  VjUi j  da  -  JtjU.  dx  (z  e  Q  a  x  e  (39) 

a  Tj 

The  essential  boundary  conditions  can  be  explicitly  included  in  the  functional  via  the  use  of 
Lagrange  multipliers;  the  functional  quantity 

Jx.(Uj  -  u.)  ds 
Tu 

can  be  subtracted  from  (39)  to  give  a  modified  functional,  fl*.  It  can  be  shown  by  variational 
calculus  that  the  Lagrange  multiplier  is  equal  to  the  boundary  traction  and  thus  the  modified 
functional  can  be  written  as 

rr(u)  =  J(|<Jijeij-ViUjjdz  -  JtUjdjL  -  Jtj(uj-u.)di 


(40) 
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Note  that  this  variational  principle  corresponds  to  the  Hu-Washizu  principle  for  the  case  where  the 
compatibility  constraint  over  the  domain  is  satisfied.  This  functional  does  not  have  an  obvious 
physical  interpretation;  however,  within  the  context  of  its  numerical  approximation  the  last  term 
might  be  thought  of  as  the  work  done  by  tractions  on  Fu  due  to  the  violation  of  essential  boundary 
conditions. 


Our  objective  is  to  express  both  functionals  in  terms  of  boundary  values.  To  attain  this  objective 
we  must  eliminate  the  domain  integration.  The  strain  energy  term  can  be  rewritten  as 


a..e..  =  a  .u. . 

>j  >j  ij  Id 

=  (0.;U.)  .  -  <J..  .U. 
U  *  d  Ud  ‘ 

=  +<i5Uj 


(41) 


Substituting  (41)  into  (39)  and  applying  the  divergence  theorem  to  the  first  term  gives 


n(u.)  « 


(42) 


The  only  approach  found  for  eliminating  the  body  force  term  is  simply  assuming  an  absence  of 
body  forces  (granted  not  a  satisfying  assumption).  With  this  assumption  the  functional  reduces  to 


n(Ui)  =  jJtiUjdi  -  jtujdi 


(43) 


Similarly  the  modified  functional  can  be  expressed  in  terms  of  boundary  values  as 


(44) 


Many  researchers  (see  e.g.  Refs  4,  7,  13, 23  and  24)  have  used  the  functional  of  equation  (43)  to 
obtain  stiffness  formulations  of  the  DBEM.  Mustoe  (Ref  14),  apparently  based  on  the  work  of 
Zienkiewicz  et  al.  (Ref  7),  used  the  modified  functional  (44)  to  obtain  a  stiffness  formulation  to  the 
IDEM  which  included  the  anificial  tractions  in  the  vector  of  unknowns.  The  modified  functional 
has  also  been  used  with  the  DBEM  when  the  interface  was  "obviously  incompatible"  (e.g.,  when 
different  "shape  functions"  were  used  for  the  two  methods).  Tullberg  and  Bolteus  (Ref  24,  1982) 
considered  yet  another  modified  functional  which  constrained  the  boundary  tractions  to  satisfy 
over-all  equilibrium. 

The  following  derivation  uses  the  boundary  form  of  the  potential  energy  functional  (43)  to  obtain  a 
stiffness  matrix  for  the  IBEM.  I  have  not  seen  this  particular  derivation  in  the  literature,  but 
literature  searches  are  inherently  incomplete.  It  actually  differs  very  little  with  the  DBEM 
derivations;  we  simply  obtain  our  nodal  displacement-traction  relation  from  the  approximation  of  a 
different  set  of  integral  equations.  The  use  of  equation  (43)  to  develop  a  stiffness  matrix  usually 
implies  the  existence  of  continuity  at  the  element  boundaries.  As  will  be  discussed  later,  we 
actually  violate  inter-element  continuity  in  the  following  derivation  but  do  so  in  good  company. 

The  analogous  derivation  *  the  DBEM  has  been  expressed  in  a  more  general  form  (see  e.g.  Ref 
23).  In  the  following  derivation  we  neither:  (1)  consider  the  essential  boundary  conditions  nor  (2) 
perform  condensation  at  the  element  level.  The  key  relation  for  the  BEM  expresses  the  nodal 
tractions  in  terms  of  the  nodal  displacements.  As  mentioned  in  the  previous  section  (assuming  no 
body  forces),  this  relationship  has  been  obtained  for  the  IBEM  by  eliminating  the  artificial  traction 
vector  from  equations  (22)  giving 

i  =  Eu  (45) 


where  E  =  FG  '  as  equation  (37). 
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Using  the  above  relationship  and  expressing  the  boundary  distribution  of  traction  and  displacement 
in  terms  of  the  so  called  "shape  functions,"  we  approximate  the  functional  of  equation  (43)  as 

n(ji)  =  Y  ]i^jN^(20  mcjO  di  E  li  -  ji^jN^(2;)  fcaLl  dx  (46) 

r  Ft 

Stadonarity  with  respect  to  each  nodal  displacement  then  gives 
Kii-  f  =  0 

where 

K  =  ^CE  +  (CE)”^) 
and 

f  =  jN^(2i)lOLldx 

An  analogous  functional  and  stiffness  matrix  derivation  exist  for  potential  problems  where  the 
vector-valued  tractions  and  displacements  arc  replaced  by  scalar-valued  sources  and  potentials, 
respectively. 

Note  that  the  matrix  in  the  quadratic  term  of  the  above  functional  is  the  K  of  the  previous  section; 
however,  its  skew  symmetric  component  contributes  nothing  to  the  value  of  the  functional  and 
stationarity  gives  symmetric  component.  The  symmetric  component  of  the  stiffness  matrix  is 
exactly  what  Brebbia  retained,  based  on  heuristic  reasoning,  from  his  direct  derivation  of  a  DBEM 
stiffness  matrix. 


(47) 

(48) 

(49) 


As  indicated  in  the  introduction,  the  symmetry  issue  has  received  significant  debate.  The  Maxwell- 
Betti  reciprocal  theorem  seems  to  imply  that  a  stiffness  matrix  should  be  symmetric  for  a  linear 
elastic  domain;  however,  this  implication  assumes  that  concentrated  nodal  forces  are  weighted 
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identically  in  obtaining  generalized  nodal  forces.  Consider  for  example  a  FEM  stiffnc  .,..  matrix 
obtained  by  a  weighted  residual  formulation  using  polynomial  shape  functions  and  trigonometric 
weighting  functions.  In  this  case  the  nonsymmetry  of  the  stiffness  matrix  does  not  imply  that  the 
reciprocal  theorem  is  violated.  For  direct  derivations  of  the  BEM  stiffness  matrix,  the  weighting 
functions  used  in  obtaining  the  generalized  nodal  forces  (29)  are  such  that  nonsymmetry  of  the 
stiffness  matrix  does  imply  violation  of  the  reciprocal  theorem  (for  concentrated  nodal  forces). 

One  could  argue  that  concentrated  nodal  forces  on  F]  have  no  meaning  since  the  original  BEM 
model  is  undefined  when  a  field  (collocation)  point  and  source  point  (concentrated  load)  coincide. 
Georgiou  (Ref  2  or  56)  has  apparently  taken  this  view  in  arguing  that  the  reciprocal  theorem  is  not 
directly  applicable  to  generalized  nodal  forces.  For  an  arbitrary  distributed  traction  symmetry  only 
implies  the  reciprocal  theorem  is  satisfied  in  some  average  sense  and  nonsymmetry  does  not 
necessarily  imply  the  reciprocal  theorem  is  violated.  However  the  argument  seems  somewhat 
sterile  if  we  think  of  a  unit  concentrated  load  as  a  uniformly  distributed  load  (with  a  resultant  of 
one)  subject  to  the  limit  of  the  distribution  area  tending  to  zero. 

From  a  practical  stand-point,  if  one  is  coupling  the  BEM  to  an  existing  FEM  program  its  likely  that 
the  architecture  of  the  existing  code  is  based  on  a  symmetric  stiffness  matrix.  Thus  symmetry  is  a 
cherished  property  and  will  be  treated  as  such  in  the  remainder  of  this  study. 

In  addition  to  lack  of  symmetry,  both  the  direct  and  variational  formulations  can  produce  stiffness 
matrices  which  violate  equilibrium  and  fail  to  allow  rigid  body  motion.  Hartman  (Ref  27,  1981) 
attributed  these  problems  to  the  approximations  made  in  obtaining  the  nodal  tracdon-displacement 
relationship  (via  E  of  equations  (37)  and  (38)  for  the  IBEM  and  DBEM,  respectively).  As 
previously  mentioned,  for  the  DBEM  the  equilibrium  problem  can  be  addressed  when  formulating 
E;  equation  (13)  is  augmented  with  two  additional  equations  and  corresponding  Lagrange 
multipliers,  which  force  the  traction  distribution  to  "identically"  satisfy  equilibrium  (Ref  13).  Most 
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approaches  to  these  problems  modify  the  stiffness  matrix  rather  than  address  the  traction- 
displacement  relationship  (Ref  24). 

Hartman  appears  (I  won't  pretend  to  have  fully  digested  his  analysis)  to  have  presented  an 
excellent  analysis  of  the  numerical  errors  which  occur  in  forming  a  stiffness  matrix  based  upon  the 
DBEM.  The  numerical  solution  of  integral  equations  can  produce  a  stiffness  matrix  with  the  above 
mentioned  undesirable  properties;  in  addition,  the  numerical  approximations  result  in  an  inherent 
incompatibility  along  Fi,  which  is  discussed  in  the  next  section. 

Discontinuity  of  Boundary  Eiement  Methods 

Regardless  of  the  formulation  there  is  a  characteristic  of  BEM  formulated  stiffness  matrices  and 
BEM  solutions  in  general  which  is  often  not  realized;  BEMs*  are  inherently  incompatible  at 
the  interface  of  homogeneous  regions. 

Derivations  of  BEMs  usually  indicate  that  the  governing  differential  equations  are  satisfied  exactly 
in  the  domain  while  the  boundary  conditions  are  only  satisfied  in  an  average  sense.  Often  this  is 
interpreted  as  meaning  the  boundary  element  "shape  functions"  constrain  the  displacement  field. 
Actually  the  meaning  has  greater  depth.  This  interpretation  is  often  a  consequence  of  invalid 
analogies  between  BEMs  and  FEMs.  The  methods  are  theoretically  related  (Refs  8,7,  and  13)  and 
certainly  the  BEM  has  borrowed  much  technology  from  the  FEM;  however,  the  accepted  use  of  the 
name  "boundary  element  methods"  for  the  numerical  solution  of  boundary  integral  equation 
methods  is  misleading.  This  interpretation  is  usually  not  enlightened  by  numerical  experience  since 
many  BEM  programs  can  not  accurately  calculate  responses  in  the  near  boundary  region. 


*  This  statement  assumes  the  BEMs  are  using  the  fundamental  singular  solution  as  presented  in  previous  sections. 
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The  "shape  functions"  used  for  the  so  called  "boundary  elements"  do  not  constrain  the 
displacement  field  but  merely  offer  an  approximation  of  the  boundary  values  to  facilitate  numerical 
integration.  That  is,  these  "shape  functions"  serve  as  a  basis  for  the  approximate  distribution  of 
boundary  values  in  the  integral  equations  --  not  for  the  boundary  response  obtained  by  these 
approximated  integral  equations.  Typically  the  integral  equations  are  satisfied  in  a  collocation 
sense  (i.e.,  at  discrete  points  along  the  boundary)  and  thus  the  boundary  conditions  are  satisfied  at 
the  collocation  points  (which  coincide  with  nodal  points  for  "continuous"  element  formulations). 
In  between  these  points,  boundary  conditions  can  be  grossly  violated.  This  is  graphically 
illustrated  in  reference  35. 

Since  the  term  "shape  function"  and  "boundary  element"  only  apply  to  the  discretization  of  the 
integral  equations  the  terms  "continuous  element"  and  "discontinuous  element"*  have  the  same 
limitation.  The  integral  equations  used  in  the  calculation  of  internal  responses  give  a  continuous 
response  throughout  the  domain;  however  for  "discontinuous  elements,"  large  response 
oscillations  can  occur  near  the  boundary  value  discontinuities.!"  The  response  with  both  elements 
is  continuous  but  their  violation  of  the  boundary  conditions  will  differ. 

Similar  to  the  displacement-based  FEMs'  satisfaction  of  equilibrium  equations  in  a  nodal  sense,  the 
BEM  (using  nodal  collocation  and  "continuous  elements")  only  satisfies  boundary  conditions  in  a 
nodal  sense.  Many  researchers  indicate  that  compatibility  between  the  BEM  and  FEM  is  easily 
satisfied  by  using  the  same  "shape  functions"  (Refs  1,  7,  13,  15,  23,  25,  26  and  58);  this  is  not 


*  "Continuous  elements"  are  characterized  by  their  C®  continuity  modeling  of  boundary  values.  "Discontinuous 
elements"  are  not  collocated  at  extreme  nodes  and  thus  can  represent  a  discontinuous  boundary  value  distribution. 

If  the  discontinuities  are  prescribed  boundary  conditions  then  these  oscillations  reflect  a  limitation  in  our  theory 
since;  1)  we  assume  displacement  continuity  and  2)  traction  discontinuity  implies  a  violation  of  equilibrium  equation 
(Ib)  (neither  the  divergence  theorem  or  mean  value  theorem  are  applicable).  Otherwise  the  discontinuities  and  the 
associated  near-boundary  errors  are  simply  a  liability  of  this  element  formulation.  The  lack  of  continuity  with 
"continuous  elements"  can  also  reduce  the  accuracy  in  the  near-boundary  region  but  usually  to  a  much  lesser  extent. 
This  was  the  motivation  for  Bilgen's  (Ref  57, 1982)  development  of  a  cubic  spline  element. 
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true.*  In  geu::^ml.  continuity  (for  a  stiffness  modeled  BEM  region)  can  not  be  satisfied  a  priori 
simply  by  the  selection  of  "shape  functions."  This  is  more  obvious  for  the  IBEM  since  the  initial 
use  of  "shape  functions"  is  for  the  approximation  of  the  artificial  traction  distribution  instead  of  the 
actual  boundary  values.  Zienkiewicz  etal.  (Ref  7)  and  Mustoe  (Ref  14)  realized  the  incompatibility 
which  existed  in  their  variational  derivation  of  an  IBEM  stiffness  matrix  based  upon  (44). 
However,  it  appears  they  attributed  this  incompatibility  to  the  use  of  the  modified  functional  since 
their  other  direct  and  variational  derivations  explicitly  indicated  that  continuity  was  satisfied.  The 
derivations  of  this  study  have  ignored  this  inter-element  continuity  violation  like  most  derivations 
in  the  literature.  This  incompatibility  is  numerically  demonstrated  for  the  IBEM  in  a  following 
section. 

» 

Incompatibility  does  not  necessarily  imply  poor  results.  On  the  contrary  many  incompatible  finite 
elements  exhibit  improved  performance,  and  many  accurate  BEM  solutions  have  been  reported. 
Recognition  of  this  basic  incompatible  behavior  could  be  of  practical  use.  It  affects  the 
convergence  properties  of  the  method  and  can  potentially  explain  other  solution  characteristics. 
One  very  useful  application  could  be  in  mesh  refinement.  Whether  the  refinement  is  a  manual  or  an 
adaptive  process,  the  variation  in  computed  responses  between  collocation  points  can  provide  a 
measure  on  which  to  base  the  refinement.  For  example,  one  could  base  the  mesh  refinement 
criteria  on  the  relative  error  between  "shape  function"  interpolated  displacements  and  integral 
equation  calculated  displacements  sampled  at  a  few  points  between  element  nodes. 

Margulies  (Ref  62,  1981)  showed  that  the  use  of  a  bounded  Green’s  function  instead  of  the  free- 
space  Green's  function  (as  used  in  this  study)  allows  the  incompatibility  to  be  eliminated.  He 


*  When  the  preliminary  report  on  this  study  (Ref  38)  was  written,  I  had  not  found  any  references  recognizing  Uiis 
incompatibility.  Among  possibly  others,  Jirousek  (Refs  59,  1978;  60;  61)  and  Margulies  (Ref  62,  1981)  had 
recognized  the  incompatibility  which  exists,  but  the  compatibility  fallacy  continues  to  propagate  throughout  the 
literature. 


indicated  that  since  we  have  freedom  in  selecting  the  geometry  of  Fj,  the  appropriate  Green's 
function  ctm  often  be  readily  obtained. 
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NUMERICAL  IMPLEMENTATION 

This  section  discusses  the  numerical  implementation  of  the  IBEM  stiffness  matrix  developed 
above.  Algorithms,  in  the  form  of  pseudo-code  outlines,  provide  detail  on  key  aspects  of  the 
implementation.  These  algorithms  have  been  implemented  in  a  research  code.  The  scope  of  this 
discussion  is  limited  to  the  calculation  of  the  stiffness  matrix;  assembly  of  element  stiffness 
matrices  is  well  documented  elsewhere  (Ref  4).  Both  the  EBEM  and  DBEM  are  addressed  below 
since  their  stiffness  matrices  (35)  can  have  the  same  form.  Rudolphi  (Ref  26)  provides  a  general 
outline  for  the  calculation  of  a  stiffness  matrix  including  the  BEM  coefficient  matrix  calculations. 
His  paper  deals  principally  with  the  potential  problem  using  the  DBEM.  Though  Rudolphi's  work 
was  not  referenced  while  implementing  the  research  code,  the  main  steps  in  the  calculations  are 
independent  of  the  application  and  the  particular  BEM. 

Due  to  the  complexity  of  BEM  and  FEM  software  systems  emphasis  must  be  placed  on  the 
modular  design  of  the  coupled  software  system.  Figure  7  illustrates  the  extent  of  the  software 
coupling  used  in  this  study.  In  the  following  discussion  I  assume  BEM  and  FEM  systems  exist 
and  require  modification.  The  stiffness  matrix  calculation  can  be  organized  as  three  tasks: 

(1)  Calculation  of  the  coefficient  matrices  (G  and  F)  and  C 

(2)  Preliminary  K  calculation,  inverse  and  matrix  product  calculations  as  shown  in 
equation  (34) 

(3)  Optional  adjustments  in  K  to  satisfy  equilibrium  and  symmetry 

The  first  task  requires  modification  of  the  BEM  program.  The  second  and  third  tasks  can  be  added 
to  the  BEM  program,  comprise  a  separate  module,  or  be  incorporated  into  an  element  routine  of  the 
FEM  program.  I  prefer  to  "weakly  couple"  the  software  systems  so  that  they  can  still  serve  as 
individual  research  and  analysis  tools.  By  this  approach  the  first  task  is  completely  performed  by 
the  modified  BEM  program.  The  second  and  third  tasks  are  performed  by  a  separate  program 
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module  which  except  for  the  three  component  matrices  of  K  (i.e.,  C,  F,  and  G)  requires  no  access 
to  the  BEM  data.  An  element  routine  is  then  added  to  the  FEM  program  which  reads  data  defining 
the  "super  element"  connectivity  and  assembles  the  BEM  K  into  the  global  system  at  the  time  of 
assembly.  Solution  of  the  FEM  system  provides  nodal  displacements  from  which  stresses  can  be 
calculated  in  both  the  FEM  and  BEM  domains. 

When  boundary  tractions  or  internal  responses  are  required  in  the  BEM  region,  the  displacements 
associated  with  BEM  nodes  must  be  retained.  For  either  BEM,  boundary  tractions  can  be  obtained 
by  the  equation:  j  =  E  ji.  For  the  DBEM  the  known  boundary  values  in  combination  with 
Somigliana's  identities  (10)  provide  internal  responses.  For  the  IBEM  the  tractions  or 
displacements  could  be  used  to  solve  a  boundary  value  problem  for  the  artificial  tractions. 
Alternatively  G'^  could  be  saved  during  the  K  calculations  giving  P  =  G‘^  £.  The  effect  of  the 
artificial  tractions  is  then  integrated  to  obtain  internal  responses  according  to  (16)  through  (18). 
The  remainder  of  this  discussion  concentrates  on  the  calculation  of  the  stiffness  matrix. 

For  smaller  problems  the  calculation  of  the  matrices  comprising  K  dominates  the  numerical  effort. 
The  coefficient  matrices  are  inherently  full  and  their  formulation  requires  extensive  numerical 
integration.  The  cost  of  their  calculation  would  appear  to  increase  with  the  square  of  the  number  of 
boundary  elements.  There  is  an  elcinent  "overhead  cost"  which  follows  this  trend;  however  for  a 
program  which  uses  a  variable  order  integration  scheme  the  cost  of  the  average  element  integration 
decreases.  This  is  due  to  the  reduction  in  the  order  of  integration  with  a  decrease  in  the  ratio  of 
element  length  to  collocation  point  distance  (Refs  35  and  36). 

Calculation  of  E  requires  an  IBEM  program  to  generate  both  F  and  G  for  a  single  homogeneous 
domain.  These  are  the  same  equations  required  to  combine  homogeneous  BEM  regions  by 
approximate  satisfaction  of  equilibrium  and  continuity  along  their  interface.  Thus  many  existing 
BEM  programs  already  have  this  capability.  The  calculation  of  two  coefficient  matrices  does  not 
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double  the  numerical  effort  because  much  of  the  overhead  in  the  numerical  integrations  is  common 
to  both  types  of  coefficients. 

In  the  DBEM  both  F  and  G  are  always  calculated.  The  only  modification  which  might  be 
required,  depending  on  the  program  design,  is  the  retention  of  all  the  coefficients.  Some 
implementations  immediately  multiply  known  boundary  values  by  corresponding  coefficients  and 
sum  the  terms  into  the  known  vector  to  form  the  system  of  equations  (15). 

Calculating  the  integration  of  the  shape  function  product  matrix,  C,  according  to  \Z6)  is  effectively 
performed  at  the  element  level.  Because  of  the  locally  based  nature  of  the  "shape  functions,"  C  is 
very  sparse.  For  quadratic  elements,  rows  associated  with  a  mid-node  have  three  nonzero  terms, 
and  those  associated  with  an  extreme  node  have  five  nonzero  terms.  In  this  implementation 
tractions  and  displacements  are  interpolated  at  the  same  order  (quadratic),  and  there  are  no 
provisions  for  traction  discontinuities.  As  a  result  M(>0  =  N(2c)  and  therefore  C  is  symmetric. 
Not  providing  for  traction  discontinuities  does  place  geometrical  constraints  on  the  BEM  region. 
However  for  BEM  regions  comprised  of  several  elements  and  few  comers,  Rudolphi  (Ref  26) 
concludes  that  traction  continuity  has  a  negligible  effect  on  the  accuracy.  For  applications  where 
the  BEM  region  is  used  as  an  "infinite  element"  the  analyst  determines  the  shape  of  the  FEM-BEM 
interface  and  thus  geometric  discontinuities  pose  no  problem. 

For  a  vector  boundary  value  problem  such  as  the  elasticity  problem,  the  two  rows  of  C 
corresponding  to  a  common  node  have  identical  terms  but  differ  in  their  position,  column.  Except 
for  the  diagonal  terms  in  C  each  term  calculated  at  the  element  level  is  complete  (i.e.,  no 
summation  is  required)  and  could  be  formally  assembled  into  four  positions  of  C.  Only  the 
diagonal  terms  of  C  which  correspond  to  shape  functions  spanning  two  elements  require  the 
addition  of  a  second  term.  Thus  it  is  very  efficient  to  both  calculate  C  and  perform  its 
multiplication  with  E  at  an  element  level.  There  is  a  maximum  of  six  unique  terms  for  a  single 
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element.  Execution  and  storage  requirements  associated  with  C  are  trivial  compared  to  F  and  G 
which  are  full  coefficient  matrices. 

The  algorithm  below  outlines  the  calculation  of  C  at  the  element  level.  The  research  code  is  written 
in  the  Modula-2  language  (Ref  63)  which  is  a  descendent  of  Pascal  and  Modula.  Those  familiar 
with  Pascal  will  recognize  the  dialect  below.  I  have  attempted  to  remove  enough  syntactical 
idiosyncrasies  from  the  code  to  permit  communication.  Calls  to  Modula-2  procedures  (analogous 
to  FORTRAN  subroutines)  are  supplemented  with  short  descriptions  of  their  purpose.  Comments 
which  do  not  replace  actual  code  are  enclosed  in  (*  *)'s.  General  file  operations  and  variable 
definitions  are  omitted.  Descriptive  variable  names  supplemented  by  comments  are  used  to  define 
variables.  Many  modem  languages  provide  for  user  defined  data  types.  Nnode  and  Mnode  are 
variables  which  can  have  the  enumerated  values  of  (a,c,b),  the  node  indices  listed  by  extreme 
nodes  first.  This  corresponds  to  a  of  Figure  4  taking  the  values  of  (1,3.2)  respectively.  These 
variables,  Nnode  and  Mnode,  serve  as  nodal  indices  to  arrays  and  promote  internal  documentation. 
More  complicated  data  types  such  as  records  have  been  converted  to  arrays  to  simplify  the 
description. 


Algorithm  1.  Shape  function  product  matrix  calculation. 


PROCEDURE  GenerateShapeProductMatrix 

(ContNode,  (*  array  of  continuous  element  node  numbers 

indexed  as  (element_number,  node__index)  *) 
numeltotal,  (*  total  number  of  continuous  elements  *) 

ShapeIntOrder)  (*  numerical  integration  order  *) 

BEGIN  (*  the  GenerateShapeProductMatrix  procedure  *) 

FOR  elnum  =  1  TO  numeltotal  DO  (*  for  each  element  *) 

(*  determine  shape  function  (N)  &  weight-Jacobian  product  values  at  each 
integration  point  *) 

FOR  Nnode  =  a  TO  b  DO  (*  an  index  on  the  node  *) 

Write  the  equation  number  for  node  ContNode [elnum, Nnode]  to  the  NMfile 
END  (*  the  Nnode  loop  *) 

FOR  intpt  =  1  TO  ShapeIntOrder  DO  {*  each  integration  point  *) 

position  =  GAUSSpt [ShapeIntOrder, intpt]  (*  obtain  the  position  of  the 
integration  point  *) 

SHAPEfunctions (position, Nvalues, intpt)  (*  calculate  the  shape  function 
values  at  "position"  and  save  in  Nvalues 
indexed  as  [intpt, node_index]  *) 

Calculate  the  Jacobian  at  "position" 

WGTdetJ[ intpt]  =  GAUSSwgt [ShapeIntOrder, intpt] *J  (*  calculate  the 
integration  weight  Jacobian  product  at  each  integration  point  *) 

END  (*  the  integration  point  loop  *) 

FOR  Nnode  =  a  TO  b  DO  (*  for  each  nonzero  shape  function  *) 

FOR  Mnode  -  Nnode  TO  b  DO  (*  for  each  unique  combination  of 
shape  function,  i.e.,  N^M  thus  symmetry  is  considered  *) 
integral  =0.0  (*  initialize  integration  *) 

FOR  intpt  =  1  TO  ShapeIntOrder  DO  (*  each  int .  pt  *) 
integral  =  integral  +  Nvalues [intpt, Nnode] * 

Nvalues [intpt, Mnode] *WGTdetJ [ intpt ] 

END  (*  the  integration  point  loop  *) 

Write  the  integral  value  to  the  NMfile 
END  (*  the  Mnode,  M  shape  function,  loop  *) 

END  (*  the  Nnode,  N  shape  function,  loop  *) 

END  (*  the  element  loop  *) 

END  (*  the  GenerateShapeProductMatrix  procedure  *) 


The  six  integration  results  are  written  to  the  NMfile  in  the  order:  a-a,  a-c,  a-b,  c-c,  c*b,  and  b-b, 
where  the  letters  indicate  the  node  index  of  each  "shape  function"  in  the  product. 

With  all  of  the  component  matrices  for  K  calculated,  the  second  task  of  calculating  an  initial 
unsymmetri-Tl  stiff  matrix  can  proceed.  The  third  task  of  satisfying  equilibrium  and  obtaining 
a  symmetric  form  can  be  combined  with  the  matrix  operations  of  the  second  task.  The  details  of 
integrating  the  two  tasks  are  highly  dependent  on  the  method  used  to  obtain  symmetry  and 
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equilibrium.  For  this  research  code,  the  tasks  are  independent,  so  different  methods  can  be 
investigated.  This  also  allows  the  symmetry  of  the  original  stiffness  matrix  to  be  examined. 

The  initial  steps  required  in  the  calculation  of  K  differ  in  the  two  BEMs.  This  difference  results 
from  the  reversal  of  F  and  G'^  in  the  traction-displacement  relations,  (37)  and  (38).  As  a  result  of 
this  difference  the  IBEM  is  burdened  with  n^  additional  multiplications  and  (n-l)n2  additions  if  E 
is  to  be  obtained  in  a  column  order.  Additionally  the  DBEM  has  reduced  memory  requirements 
since  F  can  be  read  from  disk  one  column  at  a  time.  Each  column  of  E  can  immediately  be  used  to 
calculate  the  corresponding  column  of  K.* 

As  the  order  of  the  system  of  equations  increases,  the  "calculation  of  the  inverse"  in  equation  (37) 
or  (38)  becomes  the  most  costly  step.  For  Gauss  elimination  the  cost  increases  as  n^.  A  common 
approach  (Refs  26  and  58)  to  reducing  the  cost  of  "calculating  the  inverse"  is  to  subdivide  the 
single  BEM  region  into  several  BEM  regions,  providing  many  small  stiffness  matrices  instead  of 
one  large  stiffness  matrix.  This  subdivision  introduces  additional  boundary  elements,  and  thus 
more  equations,  but  the  equations  corresponding  to  the  complete  BEM  domain  are  now  block 
banded.  Mitsui  et  al.  (Ref  58)  suggests  that  with  regard  to  efficiency  there  is  an  optimal  degree  of 
subdivision.  Subdividing  the  BEM  region  may  also  require  that  traction  discontinuities  and  the 
violation  of  equilibrium  be  treated  more  rigorously  since  the  occurrence  of  geometric 
discontinuities  and  a  relatively  small  number  of  elements  would  be  more  likely.  Kagawa, 
Yamabuchi  and  Kitagami  (Ref  64)  take  this  approach  to  the  limit;  each  boundary  element  on  Fj 
belongs  to  a  unique  homogeneous  BEM  region.  Thus  each  BEM  region  is  bounded  by  three 
boundary  elements  -  one  finite  length  element  on  Fj  and  two  "infinite  boundary  elements" 


*  The  approach  presented  here  and  used  in  the  following  algtxithm  breaks  a  cardinal  rule  of  numerical  analysis  -  an 
inverse  is  calculated.  An  alternative  exist  via  some  linear  algebra:  E  =  FG'*  <=»  EG  =  F  o  G^E^  =  F^.  In 
contrast  to  the  previous  approach,  solving  this  system  of  equations  gives  E  by  rows  instead  of  by  columns.  F  could 
be  read  from  disk  a  row  at  a  time,  but  K  would  have  to  be  kept  in-core  even  if  the  unsymmeUnc  K  was  to  be  used. 
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extending  away  from  Fj.  The  bandwidth  is  minimized  but  so  is  the  accuracy.  To  improve  the 
accuracy  in  the  region  of  principle  interest,  Q**  can  be  enlarged  to  further  isolate  Fj.  Thus  for 
comparable  accuracy  the  bandwidth  can  be  reduced  at  the  cost  of  increasing  the  number  of 
equations. 

For  both  BEMs,  a  column  of  E  is  immediately  used  to  calculate  the  corresponding  column  of  K. 
The  above  algorithm  calculates  C  at  the  element  level  and  writes  the  values  to  NMflle.  An 
algorithm  to  calculate  an  initial  K,  given  the  component  matrices,  is  shown  below.  Values  of  C 
are  read  from  NMfile,  consistent  with  the  previous  algorithm.  In  this  implementation  the  pseudo¬ 
code  given  below  is  a  separate  program  (a  Modula-2  program  module). 


Algorithm  2.  Stiffness  matrix  calculation. 


MODULE  Stif fnessMatrix 

BEGIN  (*  the  Stif fnessMatrix  module  *) 

ReadCmatrix(EqNum,C,numeltotal, order)  (*  read  C:  EqNum  contains  the 
equation  numbers  for  each  node  indexed  as  [element_number, node_index] 
and  C  contains  the  integration  of  the  shape  function  products  indexed  as 
[element_number, 0 . . 5]  *) 

ReadGmatrix  (*  read  G  *) 

ReadFmatrix  (*  read  F:  For  the  Modula-2  implementation  these  two 
procedures  cause  G  and  W  to  be  read  by  modules  involved  in  the 
calculation  of  X.  This  module  does  not  have  access  to  the  data 
of  G  and  F.  *) 

(*  Calculate  the  stiffness  matrix,  K  *) 

Initialize  K  to  zero 

FOR  col  =  1  TO  order  DO  (*  for  each  column  of  K  *) 

CalcEcol (col, Ecol)  (*  calculates  a  single  column  of  X  stored  in  Ecol  *) 
FOR  elnum  -  1  TO  numeltotal  DO  (*  for  each  element  *) 

Ccount  “0  (*  initialize  the  index  for  the  C  array  *) 

FOR  nodel  =•  a  TO  b  DO  (*  node  index  for  the  first  shape  function  *) 
roweq  *  EqNum [elnum, nodel]  (*first  row  equation  for  nodel*) 

FOR  node2  =  nodel  TO  b  DO  (*  for  each  new  combination  *) 

coleq  =•  EqNum  [elnum,  node2]  (*  first  column  equation  for  node2  *) 
Oval  =  C [elnum, Ccount]  (*  extract  the  C  term  from  the  element 
array  *) 

FOR  eqincrement  -  0  TO  1  DO  (*  for  each  nodal  dof  *) 

Crow  =  roweq+eqincrement 
Ccol  ■=  coleq+eqincrement 

(*  calc,  contributions  in  both  symmetrical  terms  *) 

K [Crow, col]  -  K [Crow, col]  +  Ecol [Ccol] *Cval 
IF  roweq#coleq  THEN  (*  not  a  diagonal  term  *) 

K [Ccol, col]  =  K [Ccol, col]  +  Ecol [Crow] *Cval 
END 

END  (*  the  eqincrement  loop  *) 
increment  Ccount 
END  (*  the  node2  loop  *) 

END  (*  the  nodel  loop  *) 

END  (*  the  elnum  loop  *) 

END  {*  the  col  loop  *) 

Adjust  K  for  equilibrium  and  symmetry 

Write  the  upper  triangular  portion  of  K  to  disk 

END  (*  the  Stif fnessMatrix  module  *) 


As  previously  noted,  due  to  symmetry  and  the  vector  nature  of  the  elasticity  problem  each  element 
C  value  (except  diagonal  terms)  has  four  positions  in  the  assembled  C.  The  "eqincrement  loop" 
above  multiplies  each  element  C  value  by  four  terms  (two  for  diagonal  terms)  in  the  column  of  E 
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and  assembles  the  products  into  the  corresponding  four  K  terms.  By  this  approach  the  sparsity  of 
C  is  exploited  and  the  matrix  is  never  assembled. 

The  last  task  in  calculating  a  stiffness  matrix  is  adjusting  K  to  satisfy  symmetry  and  equilibrium. 
The  work  of  Tullberg  and  Bolteus  (Ref  24)  provides  guidance,  based  on  numerical  studies,  for 
both  of  these  problems  .  As  previously  mentioned,  numerical  studies  have  shown  the 
unsymmetric  K  obtained  by  the  direct  approach  to  be  more  accurate  than  the  symmetrical  form 
obtained  by  a  variational  approach.  However  when  combining  the  matrix  with  an  existing  FEM 
system  a  symmetric  form  is  usually  preferable.  An  exception  can  occur  if  the  FEM  is  used  to 
model  plasticity  governed  by  nonassociative  flow  rule;  in  this  case  the  resulting  FEM  system  can 
be  unsymmetric.  In  this  study  the  stiffness  matrix  is  symmetrized  by  averaging  the  off-diagonal 
terms. 

In  the  initial  implementation  of  the  research  code,  equilibrium  considerations  have  not  been 
included.  Our  principal  interest  is  in  infinite  domain  problems  where  previously  mentioned 
techniques  are  not  applicable.  For  finite  domain  problems,  one  of  the  methods  implemented  by 
Tullberg  and  Bolteus  (Ref  24)  to  solve  the  equilibrium  problem  might  be  improved  upon.  The 
equilibrium  error  for  a  given  direction  was  corrected  by  adding  the  average  error  to  each  term. 
This  approach  did  not  consider  the  relative  magnitude  of  each  term.  A  better  approximation  might 
be  obtained  by  basing  the  error  distribution  on  the  normalized  magnitude  of  each  term  where  the 
sum  of  absolute  values  of  the  stiffness  terms  is  the  normalizing  factor. 

The  FEM  program  used  in  this  study  was  a  modified  version  of  SAPIV  (Ref  65).  An 
isoparametric  eight-node  quad  was  added  to  the  program  as  well  as  the  BEM  element  routine  for 
reading  previously  generated  stiffness  matrix  (see  Figure  7).  Since  the  BEM  region  was  not 
subdivided,  local  arrays  in  the  element  assembly  routine  had  to  be  enlarged  to  allow  assembly  of 
the  super  element. 
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NUMERICAL  RESULTS 

This  section  includes  the  numerical  results  for  tv/o  classical  problems  of  elastostatics.  Both 
consists  of  a  cavity  in  an  infinite  media  subjected  to  an  internal  pressure;  the  first  is  simple  (a 
circular  hole),  the  second  is  not  (a  crack).  For  both  problems  IDEM  analyses  were  attempted  using 
the  standard  formulation  (i.e.,  using  polynonnial  "shape  functions"  and  collocation)  and  a  stiffness 
formulation  (the  symmetric  stiffness  formulation  presented  in  previous  sections).  For  the  circular 
hole  problem  the  stiffness  formulation  gave  acceptable  results  so  a  coupled  solution  was  pursued. 

Pressurized  Circular  Hole  in  an  Infinite  Media 

The  problem  is  depicted  in  Figure  8  which  shows  the  two  boundary  clement  models  used  (a  6  and 
a  12  element  model).  Using  Airy  stress  functions  in  polar  coordinates  the  exact  solution  is  found 
to  be 


where  "a"  is  the  radius  of  the  cavity.  In  particular  I  considered  the  problem  where:  Pq  =  100,  a  = 
3,  E  =  207,900  and  v  =  0.1  (Refs  8  and  17).  It  is  interesting  to  note  that  with  this  stress  field  all 
the  strain  energy  is  distortional,  the  plane  strain  and  plane  stress  solutions  coincide,  and  the 
associated  plane-strain  problem  gives  the  exact  solution  to  the  plane  problem.  The  same  holds  true 
for  the  problem  of  a  circular  cavity  subjected  to  a  uniform  shear,  thus  these  results  hold  for  any 
uniform  load  (of  constant  direction  with  respect  to  the  polar  system).  At  the  other  extreme  we  can 
consider  the  original  problem  with  a  finite,  circular,  outer  boundary.  As  the  outer  boundary 
approaches  the  inner  boundary  (the  case  of  a  thin  walled  shell  for  plane  strain)  all  the  strain  energy 
becomes  dilational. 
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The  standard  IBEM  analysis  was  posed  as  both  a  traction  and  a  displacement  boundary  value 
problem  where  the  radial  displacement  was  obtained  from  the  a’  /e  analytical  solution.  Both 
problems  were  of  interest  because  the  stiffness  formulation  uses  the  coefficient  matrices  of  both 
problems.  The  results  for  the  six  element  model  (model  1)  are  compared  with  the  analytical 
solution  in  Figure  9.  Overall  the  results  of  the  displacement  boundary  value  problem  are  the  most 
accurate,  however  it  also  has  its  greatest  error  in  the  most  critical  region  -  near  the  boundary.  The 
results  for  the  twelve  element  model  (model  2)  are  compared  with  the  analytical  solution  in  Figure 
10.  The  accuracy  of  both  boundary  value  problems  improved  with  the  refined  discretization.  Note 
that  unlike  most  domain-based  methods  the  BEM  calculates  displacements  and  stresses  to 
comparable  accuracies. 

The  next  models  considered  for  this  problem  used  the  IBEM  in  the  previously  presented  symmetric 
stiffness  formulation.  The  IBEM  was  used  to  calculate  a  stiffness  matrix  which  was  then  read  into 
SAP  rv  to  calculate  the  nodal  displacements.  With  the  nodal  displacements,  the  IBEM  solved  the 
resulting  displacement  boundary  value  problem  for  the  nodal  artificial  tractions  and  then  any 
subsequent  internal  responses. 

The  results  for  the  six  (model  3)  and  twelve  (model  4)  elentent  models  are  given  in  Figures  1 1  and 
12,  respectively.  Note  that  the  responses  along  radial  lines  thru  element  extreme  and  mid-side 
nodes  differ.*  This  difference  indicates  the  results  oscillate  about  the  solution  with  respect  to 
change  in  the  angle.  For  the  domain  as  a  whole  there  is  significant  improvement  in  the  results  with 
the  model  refinement.  However,  while  the  displacements  improve  on  the  boundary  the  stresses  do 
not  appear  to  do  so.  This  behavior  will  be  discussed  when  considering  the  coupled  solution  that 
follows. 


*  For  Figure  1 1  response  lines  along  the  x-axis  and  y-axis  correspond  to  lines  thru  the  extreme  and  mid-side  node 
respectively. 
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Prior  to  considering  the  response  of  the  coupled  solution  let's  examine  the  degree  of  asymmetry  in 
the  stiffness  matrices  of  models  3  and  4.  In  Figure  13  a  measure  of  asymmetry  in  each  term  of  the 
stiffness  matrix  is  indicated.  A  blank  indicates  the  stiffness  value  was  less  than  10'^.  An  integer  n 
in  the  i-j  position  implies  that  KKy  -  (Kij+Kji)/2)/Kijl  <  n%  (i.e.  the  difference  in  the  old  and  new 
stiffness  terms  relative  to  the  magnitude  of  the  diagonal  stiffness  term  of  the  i***  row).  This 
measure  of  asymmetry  was  used  because  it  emphasizes  the  asymmetry  in  the  significant  terms 
(those  of  large  magnitude)  of  the  stiffness  matrix.  For  example,  an  asymmetry  of  50%  relative  to 
(Kjj+Kji)/2  is  not  very  significant  if  the  magnitude  of  Kjj  is  such  that  1000Kjj<Kji.  Due  in  part  to 
this  measure  the  asymmetries  are  "magnitude  banded"  around  the  diagonal.  The  nodes  were 
simply  numbered  in  a  counter-clockwise  manner;  so  the  large  relative  asymmetries  correspond  to 
dof  associated  with  nodes  that  are  geometrically  close.  For  both  models  the  degree  of  asymmetry 
was  on  the  order  of  5%. 

Figure  14  shows  the  results  obtained  using  a  simple  4  element  FEM  model  (model  5)  which 
exploited  the  axisymmetry  of  the  problem.  A  60  degree  sector  was  modeled  where  the  first 
element  had  a  radial  length  of  3  and  each  successive  element  doubled  in  length.  The  outer 
boundary  was  at  a  radius  of  48  and  the  nodes  on  this  boundary  were  fixed  against  translation.  The 
overall  accuracy  was  very  good.  As  expected,  note  the  improved  accuracy  of  at  the  Gauss 
points  compared  to  the  center  of  the  element.  The  numerical  results  do  not  provide  a  great  deal  of 
motivation  for  a  coupled  approach,  but  that  is  due  to  the  simplicity  of  the  problem.  Without 
axisymmetry  in  loading  and  geometry,  a  large  model  would  be  required  and  the  bandwidth  would 
necessarily  be  quite  large.  The  motivation  increases  in  dynamics  applications  where  artificial 
boundaries  can  produce  reflections  from  waves  that  should  continue  to  propagate  outward. 

The  two  coupled  models  are  shown  in  Figures  15.  For  both  models  Pj  is  a  circle  with  a  radius  of 
6.  Figures  16  and  17  give  the  numerical  results  for  the  coupled  models.  Again  the  results  indicate 
a  convergence  tendency  with  model  refinement.  As  with  the  models  3  and  4,  the  stresses  along  the 
boundary  to  the  BEM  region  do  not  appear  to  improve,  but  the  error  does  become  more  localized. 
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Even  though  the  BEM  region  has  its  maximum  error  along  Fi,  it  appears  that  this  oscillation  in  the 
response  does  not  cause  the  accuracy  in  the  FEM  region  to  be  sacrificed. 

Figures  18  and  19  summarize  the  behavior  of  the  IBEM  stiffness  formulation  along  boundaries  — 
both  by  itself  and  coupled.  As  discussed  in  a  previous  section  the  "shape  functions"  do  not  serve 
as  a  basis  for  the  boundary  response  for  the  BEM.  This  is  graphically  illustrated  in  both  figures. 
In  Figure  18,  the  response  lines  labeled  as  "BEM"  are  the  values  given  by  the  integral  equations 
rather  than  being  simply  interpolated  between  nodal  values  by  the  "shape  functions."  Due  to  the 
use  of  collocation  the  response  estimates  coincide  at  the  collocation  points  (the  nodal  points  in  this 
formulation).  The  relative  errors  at  the  extreme  and  mid-side  nodes  for  models  3  and  4  are 
Extreme  Nodes  Mid-side  Nodes 

Models  -14.3%  7.73% 

Model  4  -5.46%  3.64% 

where  a  negative  sign  indicates  an  underestimate  of  the  value.  Similarly  the  errors  for  models  6 
and  7  along  Fj  are 

Extreme  Nodes  Mid-side  Nodes 

Model  6  -5.99%  1.80% 

Model  7  -2.21%  1.21% 

The  smaller  errors  for  the  coupled  problem  are  probably  due  to  two  factors;  1)  Fj  is  in  a  region 
with  a  lower  stress  gradient  and  2)  the  FEM  region  tends  to  resist  the  distortion  of  Fj. 

For  both  the  IBEM  stiffness  formulation  and  the  coupled  approach,  the  discrepancy  between  the 
two  response  curves  for  a  given  model  gives  a  measure  of  the  goodness  of  the  discretization.  The 
"shape  function"  curve  gives  the  assumed  distribution  used  in  calculating  the  stiffness  matrix,  and 
the  BEM  curve  gives  the  resulting  response  as  calculated  by  the  integral  equation  (for  the  IBEM 
equation  [16]);  so  if  the  discrepancy  is  "significant"  the  assumed  distribution  was  meaningless. 
For  a  BEM  analysis  (standard  or  stiffness  formulation)  the  exact  solution  occurs  when  the  integral 


equation  calculations  exactly  yield  the  boundary  values.  For  a  coupled  analysis  we  can  think  of  the 
discrepancy  as  an  indication  of  the  incompatibility  between  the  two  methods.  In  this  particular  case 
the  two  subdomains  overlap  rather  than  form  a  gap. 

In  considering  the  discrepancy  between  the  assumed  distribution  of  boundary  values  and  the 
distribution  given  by  the  integral  equations,  it  is  instructive  to  consider  a  more  basic  problem. 
Consider  the  problem  of  the  circular  cavity  (geometrically  given  by  model  3)  subjected  to  a  radial 
load  of  1(XX)  at  one  node  with  the  remainder  of  the  nodes  fixed.  Figures  20  and  21  give  the  results 
to  this  problem  along  a  portion  of  F  for  the  cases  of  the  load  applied  to  both  a  mid-side  and  extreme 
node,  respectively. 

The  integral  equation  results  have  a  higher  order  of  continuity  than  the  "shape  functions,"  and  both 
the  integral  equation  curves  are  of  similar  shape.  For  the  IBEM,  "shape  functions"  are  used  to 
interpolate  the  artificial  tractions  in  the  standard  formulation  but  not  the  real  boundary  values. 
(This  differs  with  the  DBEM  where  the  unknowns  in  the  system  of  equations  are  nodal  values  of 
the  real  boundary  values  [15].)  In  forming  the  stiffness  matrix  the  IBEM  also  uses  the  "shape 
functions"  to  interpolate  the  boundary  values  (see  equations  [27],  [29]  and  [46]).  This  at  least 
partially  accounts  for  the  significant  reduction  in  accuracy  when  using  the  stiffness  formulation. 

Note  that  the  "shape  functions"  predicted  by  the  BEM  are  not  locally  based,  as  assumed,  and  are 
only  zero  at  collocation  points.  It  might  be  possible  to  use  analyses  such  as  this  (in  an  automated 
manner)  to  adjust  some  element  parameters  that  would  slightly  change  the  original  "shape 
functions"  used  to  interpolate  boundary  values.  However  it  would  probably  be  more  effective  to 
simply  use  the  discrepancy  to  determine  the  need  for  mesh  refinement.  It  is  interesting  to  note  that 
if  the  "shape  functions"  predicted  by  the  integral  equations  were  used  to  obtain  generalized  nodal 
forces  (49),  the  mid-side  node  and  extreme  node  generalized  forces  would  decrease  and  increase, 
respectively  -  exactly  what  appears  to  be  needed  (Figures  18  and  19). 
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Note  the  similarities  here  with  an  eigenfunction  problem  (especially  for  the  DBEM).  We  seek  a 
distribution  of  boundary  values  (both  tractions  and  displacements)  that  when  operated  on  by 
integral  operators  gives  back  the  same  distribution  of  boundary  values.  We  use  the  difference  in 
the  distributions  to  evaluate  the  iteration  and  continue  to  iterate  (refine  the  mesh)  if  necessary. 

For  both  the  IBEM  stiffness  formulation  and  the  coupled  solution  the  stresses  along  Fj  do  not 
improve  with  mesh  refinement.  This  is  consistent  with  the  displacement  results  shown  in  Figures 
18  and  19.  The  amplitude  and  period  of  the  oscillation  in  the  boundary  values  decrease  with  mesh 
refinement  but  in  a  manner  that  leaves  the  derivatives  of  the  displacement  (on  the  boundary)  nearly 
unchanged  relative  to  a  given  position  on  an  element.  At  practical  levels  of  discretization  this 
oscillatory  error  in  the  near  boundary  region  may  always  exist.  The  key  point  though  is  that  this 
error  had  no  apparent  effect  on  the  accuracy  in  the  FEM  region.  This  might  not  be  the  case  for  all 
problems. 

From  an  energy  point  of  view  one  might  expect  the  average  displacements  in  Figures  18  and  19  to 
under  estimate  the  theoretical  value.  Since  the  predicted  radial  displacement  is  nonuniform,  the 
model  predicts  the  presence  of  dilational  energy  as  well  as  the  expected  distortional  energy.  Since 
the  stiffness  associated  with  dilation  is  greater  than  that  of  distortion,  the  reduced  average 
displacement  is  consistent  with  this  apparent  redistribution  of  strain  energy. 

For  this  simple  problem  the  coupled  solution  performed  reasonably  well.  The  errors  are  small  and 
are  localized  around  Fj.  This  problem  suggests  that  the  coupled  solution  approach  docs  have 
potential  when  the  BEM  is  used  to  model  the  infinite  domain.  In  this  problem  we  considered  the 
best  possible  shape  for  a  cavity;  in  the  next  problem  we  consider  the  worst  possible  shape. 

Pressurized  Crack  in  an  Infinite  Media 

This  problem  consists  of  a  unit  length  pressurized  crack  in  an  infinite  media  as  shown  in  Figure 
22.  The  motivation  for  this  problem  is  the  idea  of  using  the  BEM  to  form  a  crack  tip  element  for 
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coupled  solutions.  My  particular  IDEM  stiffness  formulation  does  not  provide  for  rigid  body 
motions,  so  it  is  not  possible  to  limit  the  BEM  region  to  the  crack  tip;  instead  the  entire  crack  is 
modeled  by  the  boundary  element  method. 

The  difficulty  with  this  problem  is  that  the  stress  has  a  singularity  of  1/Vr  at  the  crack  tip.  For  this 
problem  I  chose  the  strain  energy  density  (SED)  as  a  measure  of  the  accuracy  in  the  tensor  stress 
field.  Figure  23  shows  the  analytical  solution  (based  on  Westergaard's  solution  [Ref  66])  for  the 
strain  energy  density  near  the  crack  tip  and  emphasizes  the  severity  of  the  stress  gradients. 

Fracture  mechanics  problems  are  often  solved  using  complex  series  (Ref  67).  In  the  vicinity  of  the 
crack  tip  only  the  singular  term  in  the  series  expansion  is  important.  Thus  to  evaluate  a  numerical 
solution  one  should  consider  a  region  where  the  singular  term  dominates  the  problem.  As  shown 
in  Figure  24,  the  singular  solution  is  the  major  contributor  to  the  strain  energy  density  for  the 
region  to  be  considered  in  the  subsequent  analyses. 

Figure  25  shows  the  various  geometries  of  the  IBEM  crack  models  considered  in  this  study.  To 
apply  the  IBEM  or  the  DBEM  to  a  crack  problem  requires  the  introduction  of  a  gap  (a  separation  of 
the  crack  faces)  into  the  model.  The  necessity  of  a  gap  is  easiest  to  explain  for  the  IBEM. 
Consider  a  horizontal  crack  subjected  to  mode  I  (i.e.  opening  mode)  loading.  By  symmetry,  on 
the  opposing  crack  faces  the  vertical  artificial  tractions  are  expected  to  be  of  equal  magnitude  but 
opposite  direction:  the  horizontal  artificial  tractions  are  expected  to  be  equivalent.  The  use  of  the 
term  "symmetry'  suggests  that  the  crack  faces  are  on  opposite  sides  of  the  crack  axis.  So  one 
might  consider  the  geometry  of  the  crack  as  being  defined  as  the  limit  of  two  crack  faces 
approaching  an  axis.  If  this  limiting  process  is  not  considered  when  applying  the  IBEM  the 
numerical  results  will  not  be  unique.  Consider  the  integral  equations  which  the  IBEM 
approximates  (19  and  21)  for  the  case  of  coincident  crack  faces.  Symmetry  requires  the  vertical 
tractions  on  opposing  crack  faces  to  be  equal  in  magnitude  and  opposite  in  direction.  Since  they 
are  applied  along  the  same  line  their  effects  in  the  integral  equations  cancel,  thus  their  magnitude 
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can  not  be  uniquely  determined.  Numerically  this  manifests  itself  as  a  singular  matrix  --  in  the 
reality  of  finite  precision,  a  very  ill-conditioned  matrix.  Due  to  finite  precision  the  method  gives 
results  but  they  are  absolutely  useless. 

The  introduction  of  a  gap  into  the  model  can  be  thought  of  as  our  numerical  approximation  to  the 
above  mentioned  limit.  As  shown  in  Figure  22,  for  each  model  the  elements  nearest  the  tips  taper 
to  close  the  cavity.  (Model  4  used  2  elements  on  each  face  to  taper  toward  a  tip.)  If  the  gap  is  too 
large  the  model  fails  to  geometrically  resemble  a  crack.  If  the  gap  is  too  small:  1)  the  method  is 
prone  to  roundoff  error,  2)  the  system  of  equations  can  be  very  poorly  conditioned,  and  3)  the 
numerical  integrations  can  be  come  very  expensive.  The  first  problem  occurs  with  approximately 
equal  magnitude  vertical  tractions  separated  by  a  small  distance.  At  a  response  point  relatively  far 
from  these  tractions,  the  difference  in  their  effects  (via  the  integral  equations)  would  be  many  digits 
into  the  real  word  —  thus  producing  roundoff  error.  The  third  problem  is  due  to  the  close 
proximity  of  collocation  points  to  neighboring  elements.  As  the  gap  decreases  the  collocation 
points  become  closer  to  neighboring  elements  and  thus  tlie  cost  of  numerical  integrations  increases. 
The  crack  problem  can  not  even  be  practically  attempted  with  a  standard  BEM  code  unless  it  can 
accurately  calculate  responses  in  the  near  boundary  region. 

The  following  figures  (26,  27  and  28)  summarize  a  parameter  study  that  varied  the  boundary 
discretization  (model)  and  the  gap.  The  symbols  used  for  a  given  gap  were  defined  in  Figure  25; 
in  addition,  the  length  of  the  dash  for  each  curve  indicates  the  relative  gap  sizes.  The  solid  curve 
near  the  bottom  of  each  plot  gives  the  error  in  the  strain  energy  density  for  the  singular  term 
of  the  analytical  solution  compared  to  the  complete  solution.  A  positive  error  indicates  an 
under  prediction  of  the  strain  energy  density. 

In  the  following  results  there  are  two  compounding  difficulties  for  the  IBEM  -  both  associated 
with  the  geometric  discontinuity.  The  most  obvious  is  the  infinite  stresses  at  the  crack  tip.  The 
second  problem  is  the  IBEM's  difficulty  in  dealing  with  any  geometric  discontinuity;  near 
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geometric  discontinuities  the  artificial  tractions  appear  to  become  unbounded.  Both  problems  are 
addressed  in  concentrating  elements  near  the  crack  tip. 

In  Figure  26  the  error  initially  decreases  away  from  the  tip  followed  by  a  subsequent  increase. 
Note  that  the  accuracy  is  highest  for  the  largest  gap.  At  flrst  glance  it  might  appear  that  the  increase 
in  error  is  due  to  the  previously  mentioned  roundoff  error  problem.  However  (though  not  shown) 
the  error  does  decrease  again.  Apparently  this  increase  in  error  with  distance  from  the  tip  is  a  result 
of  the  IBEM’s  difficulty  with  geometric  discontinuities.  The  initial  increase  in  accuracy  with  gap 
size  (i.e.,  as  the  geometric  discontinuity  decreases)  reinforces  this  conclusion.  For  these  results 
the  elements  adjacent  to  the  tip  are  internally  collocated  (i.e.,  a  so  called  "discontinuous" 
formulation)  to  prevent  collocation  at  the  crack  tip  where  a  traction  discontinuity  exist. 

Figures  27  and  28  further  enforce  the  above  conclusion.  In  the  region  where  I  have  attributed  the 
error  to  the  IBEM's  problem  with  discontinuities,  the  analysis  based  on  the  largest  gap  performs 
the  best.  At  a  distance  further  from  the  tip  the  smaller  gap  analyses  give  the  better  results.  Note 
that  the  smallest  gap  does  not  give  the  best  results  (see  e.g.  Figure  28).  All  of  the  previously 
mentioned  problems  which  might  occur  with  a  small  gap  are  suspect.  There  is  another  problem 
which  occurs  with  the  method  that  is  apparent  more  from  a  physical  perspective.  Recall  the  initial 
derivation  of  the  IBEM.  There  we  considered  our  actual  problem  as  embedded  in  an  infinite 
domain,  and  the  artificial  tractions  were  determined  on  the  boundary  such  that  the  boundary 
conditions  were  satisfied.  For  this  problem  embedding  the  problem  in  the  infinite  domain  implies 
that  material  exist  inside  the  crack.  As  the  gap  decreases,  the  stiffness  (relative  to  crack  opening) 
of  this  sliver  of  material  in  the  crack  increases;  that  is,  it  requires  greater  artificial  tractions  to  open 
the  crack.  For  this  problem  a  gap  in  the  range  [0.0025,0.005]  gave  the  best  results. 

Figure  29  illustrates  the  improvement  of  the  results  with  mesh  refinement.  (In  this  case  the  smaller 
the  dash  the  finer  the  mesh  near  the  crack  tip.)  As  the  mesh  is  concentrated  near  the  crack  tip  the 
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error  resulting  from  the  singularity  and  the  IBEM's  discontinuity  problem  (apparently  the  dominant 
problem)  becomes  progressively  more  localized. 

Figure  30  addresses  the  effect  of  collocating  at  the  element  crack  tip.  Tliis  analysis  was  performed 
in  preparation  for  an  IBEM  stiffness  formulation.  The  specified  boundary  values  were  zero 
traction  at  the  tip.  I  expected  a  reduction  in  accuracy  when  collocating  at  the  crack  tip;  such  was 
not  the  case.  Collocating  at  the  crack  tip  actually  gave  better  results  overall.  For  distances  greater 
than  0.03  from  the  tip,  model  3  gave  relative  errors  of  less  than  5%. 

For  this  problem,  the  need  for  a  better  discretization  was  evident  by  a  comparison  with  the  exact 
solution.  In  a  real  application  there  are  two  measures  which  could  be  considered  in  determining  the 
goodness  of  a  discretization  based  on  a  single  analysis:  1)  the  difference  between  the  "shape 
function"  interpolated  values  and  the  integral  equation  values  along  the  boundary  (as  previously 
discussed)  or  2)  an  inspection  of  the  artificial  traction  distribution.  Figures  31  give  the  artificial 
traction  distribution  along  the  upper  right  face  for  model  4  (gap  =  0.(X)5  and  "discontinuous 
elements"  are  used  adjacent  to  the  crack  tip).  As  expected  the  artificial  tractions  are  relatively 
uniform  near  the  center  of  the  crack  and  have  high  gradients  near  the  crack  tip.  The  gross  slope 
discontinuities  in  the  plot  near  *'’e  crack  tip  indicate  further  resources  must  be  concentrated  in  this 
region  to  obtain  accurate  results  near  the  boundary. 

Figure  32  gives  a  more  complete  picture  of  model  4’s  accuracy  (or  lack  thereof)  near  the  crack  tip. 
At  a  radius  of  approximately  0.05  from  the  tip,  the  error  in  the  strain  energy  density  is  15%  or  less. 
Considering  that  only  twenty  elements  were  used  to  model  the  entire  crack  and  that  the  code  was 
not  modified  to  deal  with  the  singularity,  the  results  for  this  model  are  reasonably  good. 

As  with  the  previous  problem  the  stiffness  formulation  was  less  accurate  than  the  standard  IBEM 
formulation.  Model  3  gave  the  "best  of  the  bad  results;"  for  a  response  line  from  the  crack  tip 
(perpendicular  to  the  crack)  the  error  in  the  strain  energy  density  remained  at  approximately  70% 
for  distances  in  the  range  [0.05,1.0].  The  asymmetries  in  the  stiffness  matrix  were  as  large  as 
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40%.  Both  the  factors  in  the  stiffness  matrix  calculation  (35)  are  of  questionable  accuracy.  The 
nodal  traction-displacement  relationship  (E)  is  only  accurate  if  the  artificial  tractions  are  modeled 
accurately  by  the  "shape  functions."  The  "shape  function"  product  matrix  (C)  is  only  accurate  if 
the  variation  of  the  real  boundary  values  between  nodes  is  also  accurately  represented  by  the 
"shape  functions".  Near  the  crack  tip  it  is  doubtful  that  any  of  the  boundary  values,  real  or 
artificial,  are  accurately  represented  by  the  "shape  functions." 

The  crack  problem  was  investigated  to  evaluate  the  potential  of  the  method  in  fracture  mechanics. 
Considering  the  limitations  of  the  current  implementation  familiarity  rather  than  accuracy  was  the 
main  objective.  I  have  not  investigated  the  developments  in  the  application  of  the  BEM  to  fracture 
mechanics.  However,  it  is  worth  noting  a  few  modifications  that  could  be  incorporated  in  the 
current  implementation  to  improve  the  numerical  results.  (I'll  omit  the  obvious  ones  like  increasing 
the  capacity  of  the  code  beyond  88  dof  s.) 

1)  Provide  for  rigid  body  motion  in  the  stiffness  formulation  so  BEM  crack  tip  elements  could  be 
developed.  Then  all  the  resources  of  the  BEM  could  be  concentrated  near  the  crack  tip.  To 
maximize  the  accuracy  geometric  discontinuities  should  be  accounted  for.  For  a  crack  not 
subjected  to  crack  face  loading,  many  boundary  elements  could  be  concentrated  near  the  tip  with 
the  corresponding  degrees  of  freedom  eliminated  at  the  element  level. 

2)  Incorporate  singular  "shape  functions"  at  the  crack  tip  for  the  artificial  tractions. 

3)  Consider  a  more  natural  integral  equation  formulation  for  the  crack  problem.  In  particular,  the 
displacement  discontinuity  method  is  based  on  the  analytical  solution  to  the  problem  of  a  constant 
discontinuity  in  displacement  along  a  finite  line  in  an  infinite  medium  (Ref  54).  The  formulation 
follows  the  same  development  as  the  IBEM  except  the  unknowns  are  displacement  discontinuities 
instead  of  artificial  tractions. 


For  the  reader  interested  in  a  modem  treatment  of  fracture  mechanics  by  the  BEM  the  recent  survey 
of  advances  in  the  BEM  by  Cruse  (Ref  68)  is  an  excellent  starting  point.  In  addition  Mukherjee 
(Ref  69)  dedicates  a  few  chapters  in  his  text  to  the  application  of  the  BEM  to  fracture  mechanics. 
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SUMMARY  AND  CONCLUSIONS 


Fundainental  differences  in  the  finite  element  and  boundary  element  methods  result  in 
corresponding  strengths  and  weaknesses.  Thus  a  combination  of  the  mc'.iods  using  the  strengths 
of  each  may  allow  some  classes  of  problems  to  be  solved  more  effectively.  Applications  with 
infinite  or  semi-infinite  domains  and  applications  in  fracture  mechanics  are  candidates  for  combined 
solution  approaches. 

The  direct  and  indirect  boundary  element  methods  are  the  most  common  integral  equation  methods 
in  use.  They  have  their  origins  in  classical  potential  theory  and  the  fundamental  singular  solution  is 
a  key  component  of  their  derivations.  For  both  methods  the  approximations  common  to  most 
formulations  are:  (1)  integrating  in  a  piecewise  manner,  and  (2)  solving  the  integral  equations  in  a 
collocation  sense.  The  collocation  approximation  limits  the  satisfaction  of  boundary  conditions  to 
discrete  points;  as  a  result  coupled  solutions  have  inherent  incompatibilities  at  interface  boundaries. 
While  the  "shape  functions"  serve  as  a  basis  for  the  assumed  boundary  value  disuibution,  they  are 
not  a  basis  for  the  boundary  responses.  Thus  the  prescription  of  identical  shape  functions  in 
adjoining  regions  does  not  imply  the  existence  of  compatibility. 

Coupling  approaches  can  be  categorized  by  several  parameters.  The  primary  classification  in  this 
study  indicated  which  numerical  method  hosts  the  coupling.  That  is,  is  the  final  form  of  the 
equations  "BEM  like"  or  "FEM  like"?  This  study  addresses  FEM-hosted  or  stiffness  couplings. 
These  can  be  further  categorized  by  the  derivation  of  the  stiffness  matrix.  Either  the  stiffness 
matrix  is  obtained  directly  or  by  a  variational  statement  of  the  problem.  For  many  direct  and 
variational  formulations  the  key  matrices  are  the  fully  populated  coefficient  matrices  generated  by 
the  BEMs  and  a  sparse  matrix  obtained  by  the  integration  of  "shape  function"  products.  The  BEM 
coefficient  matrices  are  manipulated  to  establish  a  relationship  between  the  boundary  tractions  and 
displacements;  this  relation  is  then  used  directly  to  obtain  a  stiffness  matrix,  or  substituted  into  a 
boundary  variational  statement. 
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For  FEM-hosted  couplings  in  elastostatics,  the  methods  can  be  combined  without  requiring  major 
modifications  to  the  existing  FEM  and  BEM  computer  programs.  This  allows  each  program  to 
serve  as  an  independent  analysis  tool  or  as  one  component  of  a  coupled  approach. 

The  major  problems  with  the  BEM  stiffness  matrices  are  their  lack  of  symmetry,  violation  of 
equilibrium,  and  inability  to  represent  rigid  body  motions.  A  variational  formulation  inherently 
produces  a  symmetric  system  which  is  equivalent  to  taking  the  symmetric  component  of  the 
corresponding  directly-formulated  stiffness  matrix.  Investigators  have  offered  several  alternatives 
for  dealing  with  the  equilibrium  and  rigid  body  problems.  They  range  from  the  introduction  of 
Lagrange  multipliers  which  enforce  equilibrium  to  ad  hoc  adjustments  of  the  stiffness  matrix. 

In  many  cases,  stiffness  matrix  calculations  for  both  BEMs  follow  the  same  main  steps.  However, 
the  direct  BEM  is  slightly  better  suited  to  the  coupled  solution  approaches  considered  in  this  study 
because  of  its  simpler  traction-displacement  relationship.  For  smaller  problems,  the  calculation  of 
matrices  comprising  the  stiffness  matrix  accounts  for  the  most  numerical  effort.  For  larger 
problems  the  "inverse  calculation"  of  a  fully-populated  matrix  dominates  the  numerical  effort.  The 
only  known  method  for  reducing  this  cost  is  to  subdivide  the  homogeneous  boundary  element 
region  (fl®)  into  a  number  of  smaller  regions— trading  the  calculation  of  a  large  inverse  for  the 
calculation  of  several  smaller  ones.  The  subdivision  of  also  reduces  the  bandwidth  of  the  final 
system  of  equations.  However,  to  maintain  comparable  accuracy  it  appears  the  interface  boundary 
(Fj)  must  be  positioned  at  a  greater  distance  from  the  region  of  interest  (usually  a  localized  region 
within  the  finite  element  region). 

Alternative  formulations  where  the  BEM  is  incorporated  by  considering  it  as  a  boundary  condition 
to  the  FEM  region  might  be  more  economical  than  the  FEM-hosted  couplings  considered  in  this 
study.  These  formulations  were  not  considered  in  detail.  I  am  not  sure  if  the  objective  of  "weak 
software  coupling"  would  be  met  by  these  formulations. 
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The  only  matrix  in  the  stiffness  calculation  not  calculated  by  standard  BEM  programs  is  the  "shape 
function"  product  matrix  C.  For  continuous  boundary  elements  C  is  symmetric.  Additionally  the 
locally  based  nature  of  the  "shape  functions"  makes  C  very  sparse.  The  calculation  of  C  and  its 
product  with  other  matrices  are  effectively  performed  at  the  element  level. 

Many  derivations  of  BEM  stiffness  matrices  have  unknowingly  ignored  the  incompatibility  along 
Tj.  Both  numerical  examples  in  this  study  indicated  that  the  standard  IBEM  is  more  Accurate  than 
the  stiffness  formulation  considered.  However,  as  shown  by  the  circular  cavity  problem 
reasonably  good  results  can  still  be  obtained  by  the  coupled  solution  except  for  a  local  region  near 
Fj  (i.e.,  the  overall  accuracy  is  not  sacrificed).  For  this  problem  the  BEM  is  used  to  model  the 
infinite  domain  -  a  role  it  is  well  suited  for.  I  am  unaware  of  any  convergence  proofs  which 
address  the  coupled  solution  approach  when  the  incompatibility  has  been  ignored.  Though 
numerous  numerical  studies  suggest  convergence  will  occur,  there’s  nothing  quite  as  comforting  as 
a  mathematical  proof.  When  the  incompatibility  is  recognized  it  can  either  be  explicitly  included  in 
a  variational  statement  (Refs  7,  59,  60,  61)  or  eliminated  by  using  a  different  Green's  function 
(Ref  62). 

The  discrepancy  between  the  assumed  boundary  value  distribution  and  that  given  by  the  integral 
equations  is  potentially  very  useful  in  mesh  refinement.  The  difference  in  the  values  could  be  used 
in  an  automated  setting  to  determine  when  to  either  subdivide  critical  elemmts  or  completely 
remesh  the  boundary. 

The  pressurized  crack  problem  served  as  the  ultimate  test  for  both  the  standard  and  stiffness 
fomiulations  of  the  IBEM.  To  address  the  problem,  the  geometrical  approximation  of  a  gap  (i.e., 
an  initial  crack  opening)  had  to  be  introduced.  In  addition  an  integration  strategy  which  is  accurate 
in  the  near  boundary  region  was  a  necessity.  While  the  accuracy  was  not  beyond  reproach,  the 
results  seemed  reasonably  good  for  a  "stock  program"  with  no  modifications  to  confront  the  crack 
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problem.  There  are  numerous  approaches  to  improving  the  accuracy,  including  completely 
different  boundary  methods  better  suited  to  the  crack  problem. 

The  inaccuracies  in  both  problems  were  interpreted  in  terms  of  the  "shape  function  fallacy."  The 
reduced  accuracy  in  the  stiffness  formulation  was  interpreted  as  being  due  to  the  increased  use  of 
the  "shape  functions"  to  interpolate  the  distribution  of  boundary  values.  The  effect  was  most 
pronounced  when  Fj  intersected  regions  of  high  gradients;  for  it  was  this  case  that  the  "shape 
function"  distribution  had  the  greatest  error.  The  only  boundary  values  we  know  are  those  at  the 
collocation  points.  This  is  brought  out  in  the  following  analogy  which  summarizes  the  results  for 
the  crack  problem. 

When  a  camper  sets  up  his  tent  he  wants  the  bottom  of  the  tent  to  remain  on  the  ground.  However, 
during  a  bad  storm  the  only  points  remaining  on  the  ground  are  those  at  the  tent  stakes.  These  are 
the  collocation  points  for  the  camper  (though  most  campers  don't  seem  to  use  this  terminology). 
When  the  "numerical  winds"  blow  the  BEM  "hangs  on"  by  its  collocation  points.  The  crack 
problem  produces  the  tornado  of  "numerical  winds."  The  tent  "beat  in  the  wind"  using  the  BEM;  it 
nearly  pulled  the  stakes  ouf  when  we  tried  to  use  the  stiffness  formulation  of  the  IBEM. 
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RECOMMENDATIONS 

The  following  recommendations  are  broken  into  two  categories:  1)  recommendations  for  future 
implementation  and  2)  recommendations  for  future  theoretical  investigations.  Implementation 
suggestions  focus  on  modifying  the  existing  IBEM  program  for  more  effective  coupled  analyses, 
and  the  development  of  other  boundary  methods  for  parallel  studies  between  the  methods.  Most  of 
the  theoretical  basis  for  this  work  is  available  in  the  literature.  Future  implementation  efforts  might 
address : 

•  an  extension  of  the  code  for  semi-infinite  domains  to  address  soil-structure  interaction 

•  symmetry  considerations  in  formulating  the  BEM  equations 

•  a  further  investigation  of  element  integrations 

-  special  techniques  for  integrations  over  singularities  (including  the  geometric 
discontinuity  case) 

-  special  techniques  for  the  near  boundary  region 

-  order  of  integration  calculations  in  variable  order  integration  schemes 

•  Galerkin  and  related  formulations  of  BEMs 

•  development  of  different  boundary  integral  method  codes  to  be  used  in  paralU  with  the 
IBEM  in  future  investigations 

-  aDBEMcode 

-  a  displacement  discontinuity  formulation  for  fracture  problems 

•  the  addition  of  infinite  boundary  elements  allowing  BEM  subdomains  to  be  subdivided 
for  coupled  solutions 

•  modification  of  the  IBEM  stiffness  formulation  to  represent  rigid  body  motions  for 
modeling  finite  domains 

•  the  addition  of  special  features  to  more  effectively  address  fracture  mechanics  problems 
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Within  the  contexts  of  elastostatics  the  following  areas  are  worthy  of  additional  theoretical 
investigation: 

•  methods  where  the  BEM  serves  as  a  boundary  condition  to  the  FEM  subdomain 

•  the  work  of  Jirousek  (Refs  59, 60,  and  61)  and  the  relationship  between  BEM 
stiffness  formulations  and  hybrid  finite  elements 

•  the  work  of  Margulies  (Ref  62)  and  alternative  Green's  functions  which  might  be 
advantageous  for  the  coupled  problems 

•  the  application  of  the  modified  functional  (40)  with  the  DBEM 

•  the  use  of  the  interface  incompatibility  in  an  adaptive  mesh  refinement  scheme 

This  last  recommendation  suggests  using  the  discrepancy  between  the  "shape  function" 
interpolated  boundary  values  and  those  predicted  by  the  integral  equations  as  a  measure  of  the 
goodness  of  a  mesh.  It  is  equally  applicable  to  BEM  and  coupled  solutions  --  "continuous"  and 
"discontinuous"  elements.  Thus  it  could  serve  as  the  needed  feed-back  for  an  adaptive  mesh 
refinement  algorithm  for  BEMs. 

Towards  our  long  term  objective  we  must  investigate  the  application  of  the  boundary  element 
method  in  elastodynamics  and  how  it  can  be  coupled  with  the  finite  element  method  in  this  context. 
Then  the  potential  of  applying  coupled  solutions  to  soil-structure  interaction  problems  can  be 
evaluated. 
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Figure  8.  Six  element  (models  1  and  3)  and  twelve  element  (models  2  and  4)  IBEM  models  of 
cuculai  cavity  problem. 
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Figure  13.  Asymmetry  in  the  stiffness  matrices  of  circular  cavity  models  3  and  4, 
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Figure  14a.  Model  5:  numerical  prediction  of  u,  versus  theory 
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Figure  14b.  Model  5:  numencal  prediction  of  o„  versus  theory 


Figure  15a.  Model  6: 6  FEs,  6  BEs,  and  SO  dof. 


Figure  ISb.  Model  7:  24  FEs,  12  BEs,  and  172  dof. 
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Figure  16a.  Model  6:  numerical  prediction  of  u,  versus  theory 


—  Theory;  m  FEM  region;  BEM  region:  x  X-axis,  +Y-axi5 

. T . ] . T . : . r . T . I 

- - - 4- - 4 - - - 4 

ww  .  t  «  I  I  I  I 


Radius 


Figure  16b.  Model  6:  numerical  prediction  of  <s„  versus  theory 
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Figure  22.  Pressurized  crack  problem. 
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Figure  25.  IBEM  crack  models. 
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28.  Model  3j  error  in  strain  energy  density  for  various  gaps. 
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Figure  29.  Error  in  strain  energy  density  for  various  models  with  a  gap  of  0.005 
(A  ~  model  1 ,  □  ~  model  2, 0  -  model  3,  x  ~  model  4). 
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Nt'R  2(1,  CO 

NF;1'.S/\  CikIc  II  K-^  (McCh'.iiic),  I’oit  Mik'Mlmik-.  CA 
NMCI)  3.  Ops  OITr;  4(1.  CO;  Ops  Dept 
NO.\.\  .loscpli  \';kIus.  Riickvlllc.  Ml) 

NORDA  CiHk-  .3(1(1,  H;i>  S(.  I.ouis.  MS 

NRl.  t  ihIc  2.s|I,  \V;isliiiii;t<iii.  DC;  CihIl'  44,'(l  ( R;inilicri:).  Washiniiloii.  DC;  Cmlc  4(i7()  (I).  r-arailas ). 

W  asliinmon,  DC 
NSC  Cotic  .S4.I,  Norlolk,  \  A 
NSD  SCr,,  Subic  Hay.  RI> 

NCSC  DC.  1  Code  '2,32  (\'ailcy).  New  l.oiuloii.  CT;  Code  44  (Carlseii).  New  l.iindoii.  Cl;  Code  44  (RS  Mimii). 

New  l.ondon,  (1;  Code  rAI,3|,  New  London.  Cl;  Lib  (Coile  45.3,3).  Newport  RI 
I’ACMISRANI  AC  111  Area.  HWO,  Kekaha,  III 

I’lllBCH  1,  CO,  San  Dieao,  CA;  1,  IkVi;.  San  Diego.  CA;  2.  CO.  Norlolk.  \A 
I’M  if  Code  KHS,  I’oim  M'lgu.  CA;  Code  .5(141.  Point  Miigu.  CA 

PWC  ACC  OITiee,  Norfolk.  VA;  Code  111.  Cireat  Lakes,  IL;  Lode  1(1.  Oakland,  CA;  Code  KH  (Library). 
Oakland,  CA;  Code  Kill.  Pearl  Harbor.  Ml;  Code  1(12.  Oakland,  C'A;  Code  12.3-0.  San  Die.no.  C  A;  Code 
.3(1.  Nortolk,  \  A;  Code  4(1(1,  Cireat  Lakes,  IL;  Code  4tltl.  Oaklani.1.  C.A;  Code  4(1(1.  P.'arl  Harbor.  HI;  Code 
40(1.  San  Diego,  C.A;  Code  420.  Cireat  Lakes.  II  ;  Coile  420.  Oakland.  C.A;  Code  42013  (Waid).  Snbie  Hay. 
RP;  Code  421  (Kayal,  Pearl  Harbor.  HI;  Code  421  (Ouin).  San  Die.no.  C'A;  Code  421  (Reynolds).  Sin 
Diego,  CA;  C  ode  422.  San  Dieno.  CA;  C  ode  42,3.  San  Diego.  CA;  Code  424.  Norlolk.  \  A;  C  ode  .5(1(1.  Cireat 
Lakes,  IL;  Code  500.  Oaklaiul.  C.A;  Library  (Code  1.34).  Pearl  Harbor.  HI;  Library,  (iiiani.  Mariana 
IslaiKb.  Library.  Norlolk,  \.A;  Library.  Pensaeola.  FL;  Library,  ^'okosuka.  lapan;  Teeli  Librtirv.  Subie  Hay. 
RP 

Spec  PWO  (Code  0S\).  Meehaiiiesburg,  P.A 

SCH.ASF.  Hannor.  PWO  (C'<H(e  S.32.3).  Hrenierton.  WA 

SL  PSlllP  leeh  Library.  Newport  News.  \'A 

OS  Dl.l-r  OF  INILRIOR  Natl  Park  Sve,  RMR  PC.  Denver.  CO 

I'SD.X  F.xt  Ser\  (1  Maher).  Washington.  DC;  For  Sve  Reg.  S.  (Bowers).  .Ail.inta,  OA;  For  Sve.  Reg  Bridge 
Lngr,  Aloha,  OR,  For  Sve.  leeh  Fngrs.  Washington.  DC 
eSNA  Ch,  Meeh  Fngrg  Dept.  .Annapolis,  MD;  Oeean  F.ngrg  Dept  ( .MeC  oriniek ).  .\nn;ipolis.  MD;  PWO. 
Aiintipolis.  MD 

CAL  POl.A'  Si  A  IT  (.(NIVFRSI  lA'  CT;  Dept  (Penalba).  San  Luis  Obispo.  CA 
CALIFORNIA  S  FA  I  F,  LNI\'F.RSi  rY  C.V.  Chelapali.  I.ong  Heaeh,  C  A 
CArilOl.lC  CNI\'  of  Am.  CF:  Dept  (Kinil.  Washington,  DC 
CUT'  OF  HFRKFLFA'  I’W.  Lngr  Dis  I  Harrison),  Berkeley.  CA 
C  i  rV  OF  l.IN  FRMORF  Daekins,  PF.  Livermore.  C  A 
CLARKSON  (01.1.  OF  TFCH  CF;  Dept  (Batson),  Potsdam,  NY 
COLORADO  S  FA  I  F  CNI\ TRSi  r'l  C  F.  Dept  (Criswell).  Ft  Collins.  ( () 

CORNF.l.l.  CNIS'F.RSIT)'  Civil  A;  Fnviron  Fngrg  (Dr.  Kulhawy).  Ithaea.  NS  ;  Library,  Ithaea.  NS' 

DAMFS  sS;  MOORF  Library,  Los  Angeles.  CA 

FLORIDA  All. .ANTIC  TiNIS’FRSTTS'  Oeean  Fngrg  Dept  (Su).  Boea  R;iton.  1-L 
FLORIDA  INST  OF  1  FC  H  C  F  Dept  (Kalajian).  Melbourne.  FL 

CiFORClA  INS'FFFfLFF  OF'  TFCHNOLOCiS'  C  F  Seol  (Kahn).  Atlanta.  CiA;  CF  Seol  (Swanger).  Atlantii.  (i.A; 

CF.  Seol  (/.uruek).  Atlant;i,  CiA.  Meeh  Fngrg  (Fulton).  Atlanta.  (i.A 
INSnrCTE  OF  MARINE  SCIFNCFS  Library.  Port  Aransas.  IX 
lOHNS  HOPKINS  l.!NIV  CF  Dept  (.tones),  Baltimore.  MD 

l.AWRENCF  LIVFRMORF.  NATL  LAB  F,l  Fokar/.  Livermore.  C.A.  Pkint  Fngrg  Lib  (l.-iod),  Livermore.  CA 
LFHIOH  TtNIVFRSI  I  S'  Linderm;m  Library.  Bethlehem.  P.A 
LONCi  BFACH  PORI  F.ngrg  Dir  (.Allen).  Long  Beaeh.  CA 
MlCHlCiAN  TFCH  CNIS'FRSTTS  CF  Dept  (Haas),  Houghton.  Ml 

MTT  F.ngrg  Lib.  Cambridge.  MA;  Lib,  Teeh  Reports,  C  ambridge.  M.A;  RS'  '('hitman.  (  ambiidge,  MA 
NATL  ACSADFMY  OF  .SX'IFN'C  FS  NR(  .  Naval  Studies  Bd.  Washington,  D( 

NFW  MEXICO  SOLAR  FNFRCiS'  INST  Dr  /.wibel.  Las  (  ruees,  NM 

NEW  YORK-NFW  IFRSFY  PORT  AfCTH  RiXD  Fngr  (S'onlar),  ,lerse\  Citv.  N,l 

NORTIIWF.STFRN  TUNIS'  C  F  Dept  (Belytsehko).  Fvaiiston,  IL 

OHIO  ST.ATF  I 'NIYFRSTTS'  CF  Dept  (H.  Aileli).  Columbus,  OH;  CF  Dept  (Sier;ikowski).  Colunibus,  OH 
ORFdON  ST.VTF  CNIS'FRSTTS'  (  F.  Dept  (Hieks).  C  orvallis.  OR.  CF  Dept  (Hudspeth),  Corvallis,  OR;  Cl 
Dept  (LeonartI).  Corvallis.  OR 

PFNNSS'I.S'.ANIA  STATE  CNIS'FRSTTS'  Areh  Fngrg  Dept  (Ciesehwiiulner),  Cniwrsity  Park,  P.A;  (iotolski. 

liniversity  Park.  PA,  Rseh  Lab  (Snyiler).  Stiile  College,  PA 
POR  TLAND  S'TA  I  F  CNIS'FRSTTS’  F.ngrg  Dept  (Miglionl.  Portlaiul.  OR 
IHBTDUF  L^NIS'FRSTTS'  CF  Seol  (Leonards),  W.  l.alayelle.  IN;  Fngrg  Lib.  SS'.  l.alayette.  IN 
RCTCiFRS  CNIS'FRSTTS'  CF,  Dept  (llanaor).  Piseataway,  JH 
S.AN  DlFCiO  PORT  Port  Fae.  Proj  Fngr.  San  Diego.  CA 
S.AN  DIFCiO  S'TA  !  T;  UNIS'  CF  Dept  (Krishnamoorthy).  San  Diego.  (  A 
SFAI'TLF  CNTv  TP.STTY  C  Fi  Dept  (Sehwaegler),  Seattle.  SVA 
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SOUTHERN  ILLINOIS  UNIVERSITY  CE  &  Mech  Dept  (Kassimali).  Carbondalc,  IL 

SOUTHWEST  RSCH  INST  Energetic  Sys  Dept  (Esparza).  San  Antonio.  TX:  King.  San  Antonio.  TX;  M. 

Polcyn.  San  Antonio,  TX;  Marehand.  San  Antonio,  TX 
STANFORD  UNIVERSITY  App  Mech  Div  (Hughes).  Stanford,  CA 

STATE  UNIVERSITY  OF  NEW  YORK  CE  Dept  (Reinhorn)  Buffalo.  NY;  CE  Dept,  Buffalo.  NY 
1EXAS  A&l  UNIVERSHY  Civil  &  Mech  Engr  Dept.  Kingsville.  TX 

TEXAS  A&M  UNIVERSITY  CE  Dept  (Machemehl),  College  Station.  TX;  CE  Dept  (Nicdzwceki).  College 
Station.  I'X;  Ocean  Engr  Proj.  College  Station.  TX 

UNIVERSHY  OF  CALIFORNIA  CE  Dept  (Fenves).  Berkeley.  CA;  CE  Dept  (Gerwick),  Berkeley.  CA;  CE 
Dept  (Herrmann).  Davis,  CA;  CE  Dept  (Kutter),  Davis,  CA;  CE  Dept  (Romstad),  Davis.  CA;  CE  Dept 
(Shen).  Davis,  CA;  CE  Dept  (Taylor).  Berkeley.  CA;  CE  Dept  (Taylor).  Davis,  CA;  Engrg  (Williamson), 
Berkeley.  CA;  Geotech  Model  Cen  (Cheney),  Davis.  CA;  Mech  Engrg  Dept  (Armand).  Santa  Barbara.  CA; 
Mech  Engrg  Dept  (Bayo),  Santa  Barbara,  CA;  Mech  Engrg  Dept  (Bruch).  Santa  Barbara.  CA;  Mech  Engrg 
Dept  (Goldsmith).  Berkeley,  CA;  Mech  Engrg  Dept  (Mitchell).  Santa  Barbara.  CA;  Mech  Engrg  Dept 
(Tulin).  Santa  Barbara.  CA;  Naval  Arch  Dept,  Berkeley.  CA 
UNIVERSITY  OF  COLORADO  CE  Dept  (Hon-Yim  Ko).  Boulder.  CO 
UNIVERSITY  OF  CONNECTICUT  CE  Dept  (Murtha-Smith),  Storrs.  CT 

UNIV'ERSITY  OF  FLORIDA  CE  Dept  (Townsend),  Gainesville.  FL;  Engrg  Sei  Dept  (Malvern).  Gainesville, 

FL 

UNIVERSITY  OF  HAWAII  CE  Dept  (Chiu),  Honolulu.  HI;  Manoa.  Library.  Honolulu.  HI;  Ocean  Engrg 
Dept  (Ertekin),  Honolulu.  HI 

UNIVERSITY  OF  HOUSTON  CE  Dept  (K.J.  Han).  Houston.  TX 

UNIVERSITY  OF  ILLINOIS  CE  Lab  (Abrams),  Urbana.  IL.  CE  Lab  (Pecknold),  Urbana.  IL;  Library. 

Urbana,  IL;  M.'T.  Davisson.  Urbana,  IL;  Metz  Ref  Rm.  Urbana,  IL 
UNIVERSITY  OF  MARYLAND  CE  Dept  (Goodings).  College  Park,  MD;  CE  Dept  (Wolde-Tinsae).  College 
Park.  MD 

UNIVERSITY  OF  MICHIGAN  CE  Dept  (Richart),  Ann  Arbor.  MI 
UNIVERSITY  OF  NEBRASKA  Polar  Ice  Coring  Office.  Lincoln.  NE 

UNIVERSITY  OF  NEW  MEXICO  HL  Schreycr,  Albuquerque.  NM;  NMERI  (Bean).  Albuquerque.  NM; 

NMERI  (Falk).  Albuquerque,  NM;  NMERI  (Leigh).  Albuquerque.  NM 
UNIVERSITY  OF  PENNSYLVANIA  Dept  of  Arch  (P.  McCleary).  Philadelphia.  PA 
UNIVERSITY  OF  SOUTH  C:AR0L1NA  Engrg  Coll  (Karabalis).  Columbia.  SC 

U'NIVERSITY  OF  TEXAS  CE  Dept  (Stokoe),  Austin,  TX;  CE  Dept  (Thompson).  Austin,  TX;  Construction 
Industrv  Inst.  Austin,  TX;  ECJ  4.8  (Breen),  Austin,  TX 
UNIVERSITY  OF  WASHINGTON  CE  Dept  (Mattock).  Seattle.  WA 

UNIVERSITY'  OF  WISCONSIN  CE  Dept  (Peyrot).  Madison,  Wl;  Great  Lakes  Studies  Cen.  Milwaukee.  W1 

WASHINGTON  DHHS.  OFE/PHS  (Ishihara).  Seattle.  WA 

/MMNA  ENGRG,  INC  Walczak,  Watertown,  MA 

ADVANCED  COMPOSITE  ENCiRG  Leewood,  West  Lafayette.  IN 

ADVANCED  TECHNOLOGY,  INC  Ops  Cen  Mgr  (Bodnar).  Camarillo.  CA 

AMERICAN  CONCRETE  INSTITUTE  Library,  Detroit,  MI 

AMETEK  OFFSHORE  RSCH  Santa  Barbara,  CA 

APPLIED  RSCH  ASSOC.  INC  Higgins,  Albuquerque.  NM 

ARM.STRONG  AERO  MED  RSCH  LAB  Ovenshirc.  Wright-Patterson  AFB.  OH 

ARVID  GRANT  &  A.SSOC  Olympia,  WA 

ASSOCIATED  SCIENTISTS  Dr.  McCoy.  Woods  Hole,  MA 

ATLANTIC  RICHFIELD  CO  RE  .Smith,  Dallas.  TX 

HATTELLE  D  Frink.  Columbus,  OH 

BECHTEL  CIVIL,  INC  Woolston,  San  Francisco.  CA 

BETHLEHEM  STEEL  f.'O  Engrg  Dept  (Dismuke).  Bethlehem.  PA 

BRITISH  EMBASSY  Sei  &  Tech  Dept  (Wilkins).  Wa.hington.  DC 

BROWN  &  ROOT  Ward.  Houston,  TX 

CANAIY.A  Uiokinetics  &  Assoc.  LTD,  Ottawa.  Ontario;  Concordia.  CE  Dept  (Marsh),  Montreal.  Oucbec; 

Viateur  De  Champlain.  D.S.A  ,  Mataite.  Tainida 
C  HEVRON  OIL  FLD  RSCH  CO  Strickland,  La  Habra.  CA 
(  llll.DS  ENGRG  C  CtRP  K.M.  Childs.  Jr,  Medlield.  MA 
CLARENCE  R  JONES.  CONSLT  IN.  LID  Augusta.  CiA 
COLLINS  ENGRCi,  INC  M  Ciarlich.  Chicago,  IL 
CONRAD  ASSOC  Luisoni.  \  ;in  Nuys.  CA 
CONSOER  TOWNSEND  sk  ASSOC  Schramm,  Chicago.  IL 
CONSTRUCTION  TEC  H  LABS.  INC  fi.  Corley.  Skokie.  IL 
1)11.1. INGHAM  C'ONSTR  C ORP  (HDikC).  F  McTIale.  Honolulu.  HI 
DKAV'C)  CORP  Wright.  Pittsburg.  PA 

TA  ALUATION  assoc.  INC'  MA  Fedele.  King  of  '-tussia.  PA 
GEIGER  ■  KKBNA  K.P.  Hampton.  Bellingham  'tVA 
Cil  N  MOTORS  RSCH  LABS  Khalil,  Warren.  Ml 


105 


(iRL.’MMAN  Af-ROSPAl  i;  CORP  IVdi  Info  C  lr.  »cthpai:c.  \V 

IIAVNI.S  cV  ASSOC  II.  II, ones.  PI-.  Oaklaiul,  CA 

IlIRSCII  eC  CO  I.  Ilirsch,  San  nicao.  CA 

11.1  DlXil-NKOl.H  ASSfH  \V  ,\Iuiiloni:h.  San  I  .ancisa).  CA 

111  (illliS  AIR(  RAI  1  CO  IVdi  Hoc  Con.  I'.l  Sou -mlo.  CA 

IIIAVAI  01  PNORS  Mil  IliiwaMt.  Oonvor.  CO 

II  I  LAND  KAVANACill  \VA  I  liRlfCRS  .  P(  Nov.  York.  NY 

IN  11.  MARIIIMP.  INC  0  Walsh.  San  Podro,  (  A 

IRl.-rnO  liipnl  Proo  Oir  (R.  Oanlonll.  laiaan.  MN 

lOIIN  ,1  MC  MCl.l.l-N  ASSfK  I.iliran.  Now  York.  NY 

l.i:o  .\  OAl.Y  CO  llonolnlu.  Ill 

M.Y  /I.  I  LIN  ASSOC.  INC  DA  Cnooo.  Now  York.  NY 

I.IN  OI  I  SIIORP.  LNCIRCI  P.  Chow.  San  I  rancisoo  CA 

LINDA  IlAl  1.  LIHRARY  Doo  Dope  Kansas  (  iii .  MO 

l.OCKlIia.D  Rsoh  l.ah  (Nonr-Ornirl).  Palo  Alto.  CA 

l.O(  KWOOD  CiRLLNi:  PNORS,  INC  KV  Ifomlaptuli,  Dallas.  I\ 

MARALllON  OH.  CO  (iamhlo.  Houston.  1  .\ 

.MARC  ANALYSIS  RSC  1 1  C ORP  Hsu.  Palo  Alio.  C  A 
M.ARMT.CH  I  NORO  Donoahuo.  .Austin.  LX 
MC  Cl.Ll.LAND  liNCiRS,  INC  l.ihiaiv.  ILuislon,  l\ 

MOlilL  R  sk  1)  CORP  Ollshoio  Lnal.a  Lih.  Dall.is.  1\ 
liDWARD  K  NODA  A;  ASSOC  Honolulu.  HI 
NLAV  /li.AL.A.ND  N/  Conorolo  Rsoh  Assoo.  l.ibrars.  Poriina 
NCllN  A;  ASSOC  AC  Nuhn.  Way/ala.  NM 
PACIFIC  MARINI.  ILC  11  (M.  VVaanorl  Duvall.  WA 
PlLl:  IR.'CK.  INC  Sinool,  Jupilor.  FI. 

PMB  SA'S  ENORO,  INC  Boa.  San  Franoisoo.  CA:  S.  Nour-OntkI.  San  F'ranoisoo.  CA 

PRESNEl.l.  ASSOC  ,  INC  DO  Prosnoll.  .Ir.  1  ouisvillo.  KY 

SANDIA  LABS  Librarv,  Livonnoro.  CA 

SARCiENT  sV  HI  RKES,  INC  .IP  Pioroo,  .Ir.  Now  Orloans.  LA 

SAL  Dl  ARABI  A  Kintj  Sauol  Hniv.  Rsoh  Con,  Riyarih 

SEA  FECI  1  CORP  Poioni,  Miami,  11. 

SHIiLL  OIL  CO  I.  Doylo,  Housion.  IX 

SIMPSON.  OL'MPliRIZ  ok  HEOER,  INC  E  Hill.  CE.  Arlington.  MA 

SRI  IN'LL  Fmarg  Mooli  Dopt  (Ciraitl),  Monlo  Park.  CA;  Engrg  Mooli  Dopi  (Simons).  Monlo  Park.  CA;  J.L. 

.lonos.  Chom  L.nai  Lab.  Monio  Park.  CA 
Llir:  1  ANLO//.1  CO.  INC  W.  Fanto/i'i.  Cuporlino.  CA 

TRW  INC  Crawlonl.  Rorlonolo  Boaoh.  CA;  F)ai.  San  BornarrIino.  CA;  Enai  Librars.  ClovolaiHl.  OH;  M  Katona, 
S;m  BornarrIino,  C  A;  Rorigois,  Rorlonrio  Boach  CA 
ICDOR  ENORO  CO  Ellogoorl.  Phooni.'r,  AZ. 

L'NLLEl)  KINODOM  I'niv  Coll  Swamsoa  (Zionkiowic/.).  Wales 

\'SE  C)oo;in  Enarg  Op  (Murton),  Aloxanrlria,  VA 

WEIDLINOER  ASSOC  F.S.  Wong,  Palo  Alto.  CA 

WL.l.LSPRlNO  COMM  11  Zarooor,  Marshall,  VA 

WES  I  INOHOCSi;  ELECLRIC  CORP  Librarv,  Pittsburg.  PA 

WISS.  .I.ANNEY.  LiLS'INER.  rk  ASSOC  DW  Pfoilor,  Northbrook.  11, 

WOODWARD-C  l  A  DE  CONSLU.  I  AN  1  S  R  Dominguo/.  Houston.  LX.  WosI  Rog.  Lib,  OaklanrL  CA 
BIRNSLIEL,  C  lorrosl  Hills,  NY 
BROWN.  ROBl  RL  I  nivorsity.  Al. 

HI  LLOCK,  11.  La  C  anaria.  C  A 
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INSTRUCTIONS 


The  Naval  Civil  Engineering  Laboratory  has  revised  its  primary  distribution  lists.  The  bottom  of  the 
label  on  the  reverse  side  has  several  numbers  listed.  These  numbers  correspond  to  numbers  assigned  to 
the  list  of  Subject  Categories.  Numbers  on  the  label  corrosponding  to  those  on  the  list  indicate  the 
subject  category  and  type  of  documents  you  are  presently  receiving.  If  you  are  satisfied,  throw  this  card 
away  (or  file  it  for  later  reference) 

If  you  want  to  change  what  you  are  presently  receiving: 

•  Delete  -  mark  off  number  on  bottom  of  label. 

•  Add  -  circle  number  on  list. 

•  Remove  my  name  from  all  your  lists  -  check  box  on  list. 

•  Change  my  address  -  line  out  incorrect  line  and  write  in  correction  (DO  NOT  REMOVE  LABEL). 

•  Number  of  copies  should  be  entered  after  the  title  of  the  subject  categories  you  select. 

Fold  on  line  below  and  drop  in  the  mail. 

Note:  Numbers  on  label  but  not  listed  on  questionnaire  are  for  NOEL  use  only,  please  Ignore  them. 


Naval  Civil  Engineering  Laboratory 
Port  Hueneme.  CA  93043-5003 


Official  Business 

Penalty  for  Private  Use,  $300 


BUSINESS  REPLY  CARD 

FIRST  CLASS  PERMIT  NO.  69 

POSTAGE  WILL  BE  PAID  BY  ADDRESSE 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


Commanding  Officer 
Code  L34 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme.  California  93043-5003 


DISTRIBUTION  QUESTIONNAIRE 

The  Naval  Civil  Engineering  Laboratory  Is  revising  its  Primary  distribution  lists. 


SUBJECT  CATEGORIES 

1  SHORE  FACILITIES 

2  Construction  methods  and  materials  (Including  corrosion 

control,  coatings) 

3  Waterfront  structures  (maintenance/deterioration  control) 

4  Utilities  (Including  power  conditioning) 

5  Explosives  safety 

6  Aviation  Engineering  Test  Facilities 

7  Fire  prevention  and  control 

8  Antenna  technology 

9  Structural  analysis  and  design  (including  numerical  and 

computer  techniques) 

10  Protective  construction  (including  hardened  shelters. 

shock  and  vibration  studies) 

11  Soil/rock  mechanics 

14  Airfields  and  pavements 

15  ADVANCED  BASE  AND  AMPHIBIOUS  FACILITIES 

16  Base  facilities  (Including  shelters,  power  generation,  water 

supplies) 

17  Expedient  roads/airfields/bridges 

18  Amphibious  operations  (including  breaKvaters.  wave  forces) 

19  Over-the-3each  operations  (Including  containerization. 

materiel  transfer,  lighterage  and  cranes) 

20  POL  storage,  transfer  and  distribution 


TYPES  OF  DOCUMENTS 

85  Techdata  Sheets  86  Technical  Reports  and  Technical  Notes 
83  Table  of  Contents  &  Index  to  TDS 


28  ENERQY/POWER  GENERATION 

29  Thermal  co.nservatlon  (thermal  engineering  of  buildings.  HVAC 

systems,  energy  loss  measurement,  power  generation) 

30  Controls  and  electrical  conservation  (electrical  systems. 

energy  monitoring  and  control  systems) 

31  Fuel  flexibility  (liquid  fuels,  coal  utilization,  energy 

from  solid  waste) 

32  Alternate  energy  source  (ge>  hermal  power,  photovoltaic 

power  systems,  solar  systems,  wind  systems,  energy  storage 
systems) 

33  Site  data  and  systems  Integration  (energy  resource  data. 

energy  consumption  data.  Integrating  energy  systems) 

34  ENVIRONMENTAL  PROTECTION 

35  Solid  waste  management 

36  Hazardous/toxic  materials  management 

37  Waste  water  management  and  sanitary  engineering 
36  Oil  pollution  removal  and  recovery 

39  Air  ^'ollution 

44  OCEAN  ENGINEERING 

45  Seafloor  soils  and  foundations 

46  Seafloor  construction  systems  and  operations  (including 

diver  and  manipulator  tools) 

47  Undersea  structures  and  materials 
46  Anchors  and  moorings 

49  Undersea  power  systems,  electromechanical  cables. 

and  connectors 

50  Pressure  vessel  facilities 

51  Physical  environment  (Including  site  surveying) 

52  Ocean- based  concrete  structures 
54  Undersea  cable  dynamics 

82  NCEL  Quides&  Abstracts  [n  None- 
91  Physical  Security  ^  remove  my  name 


